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EXECUTIVE SUMMARY 


LSPRAY-III is a Lagrangian spray solver devel- 
oped for application with parallel computing and un- 
structured grids. It is designed to be massively paral- 
lel and could easily be coupled with any existing gas- 
phase flow and/or Monte Carlo Probability Density 
Function (PDF) solvers. The solver accommodates 
the use of an unstructured mesh with mixed elements 
of either triangular, quadrilateral, and/or tetrahedral 
type for the gas flow grid representation. It is mainly 
designed to predict the flow, thermal and transport 
properties of a rapidly vaporizing spray because of 
its importance in aerospace application. The manual 
provides the user with an understanding of various 
models involved in the spray formulation, its code 
structure and solution algorithm, and various other 
issues related to parallelization and its coupling with 
other solvers. With the development of LSPRAY-III, 
we have advanced the state-of-the-art in spray com- 
putations in several important ways. Some of the 
major features of the spray module are: 

• It facilitates the use of both unstructured grids 
and parallel computing and, thereby, facilitates 
large-scale combustor computations involving 
complex geometrical configurations. 

• In order to deal with modern gas-turbine fuels 
that are mixtures of many compounds, we have 
extended the spray formulation to the modeling 
of multicomponent liquid fuels. 

• In order to extend the applicability of the spray 
computations over a wide range of low-pressure 
conditions, we have completed the implemen- 
tation of the liquid-phase variable thermo- 
transport properties into our spray formulation. 

• In order to reduce some uncertainty associated 
with the specification of the initial droplet con- 
ditions, we have completed the integration of an 


atomization module into our spray solution pro- 
cedure. The atomization module contains the 
following primary atomization models: (a) sheet 
breakup, (b) air blast, (c) blob jet, and (d) BLS 
(Boundary-Layer Stripping), and the following 
secondary droplet breakup models: (a) Rayleigh- 
Taylor, (b) TAB (Taylor Analogy Breakup), and 
(c) ETAB (Enhanced Taylor Analogy Breakup). 

• Because of its importance in some NASA-related 
critical applications, we have recently com- 
pleted the integration of a superheat vaporiza- 
tion model into our spray solution procedure. 

• The spray module is designed in such a way so 
that it could easily be coupled with other CFD 
codes. 

• It can be used in the calculation of both steady 
as well as unsteady computations. 

• We have developed and implemented a time- 
averaging method into the calculation of the 
liquid-phase source contributions in order to ac- 
celerate convergence towards a steady state so- 
lution. 

• The spray module has a multi-liquid and multi- 
injection capability. 

With the aim of improving the overall solution 
procedure involving sprays, we have made several rel- 
evant contributions to the gas side of the computa- 
tions: 

• In order to demonstrate the importance of chem- 
istry/turbulence interactions in the modeling 
of reacting sprays, we have extended the joint 
scalar Monte Carlo PDF (Probability Density 
Function) approach to the modeling of spray 
flames, unstructured grids, and parallel comput- 
ing. 
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• In order to account for the nonideal gas behav- 
ior under critical and supercritical conditions, 
we have completed the integration of a high- 
pressure EOS (Equation-of-State) into our CFD 
module. Also, we have completed the implemen- 
tation of a high-pressure correction into the cal- 
culation of the gas-phase transport properties. 

The spray solution procedure provided favor- 
able results when applied to the modeling of sev- 
eral different spray/gaseous flames - representa- 
tive of those encountered in gas-turbine combus- 
tors, stratified-charge rotary combustion (Wankel) 
engines, supersonic diffusion flames, and pulse det- 


onation combustion devices. Its use has been demon- 
strated in various important NASA projects: Ultra- 
Efficient Engine Technology (UEET), Pulse Deto- 
nation Combustion Technology (PDCT), & Rotary 
Combustion Engine Technology Enablement Project 
(RCETEP). 
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NOMENCLATURE 


a 

liquid jet radius, m, or 

m 

rh k 

rilk, flash 

CL 0 

parent drop radius, m 
initial droplet radius, m 

rilko 


outward area normal vector 


of the nth surface, m 2 

riik,t 

A 

a constant 

Bk 

Spalding mass transfer number 

m 0 

B t 

Spalding heat transfer number 

m(t) 

B 0 

a constant 

M 

Cfreq 

a constant 

M a 

c 

a constant 


c D 

drag coefficient 

M f 

Mi 

c p 

specific heat, J / (kg K) 

D 

turbulent diffusion coefficient, m 2 /s 

Dab 

diffusion coefficient, cm 2 /s 

M sbe d 

d 

droplet diameter, m 

MMD 

db 

droplet diameter after primary breakup, m 

Tlk 

Al 

ligament diameter, m 

Tl parent 
product 

n(t) 

do 

orifice diameter, m 

dparent. 

parent drop diameter, m 

dproduct 

product drop diameter, m 

N 

dt 

time increment, s 

N n 

d6 

half-cone angle, deg. 

o 

Nt 

Dparent 

energy contained in the parent 

J 


drop, (kg m 2 )/s 2 

N v 

Dproduct 

energy contained in the product 

P 

Nu 


drop, (kg m 2 )/s 2 

P 

h 

specific enthalpy, J/kg, 

Pc 


liquid sheet thickness, m 

Pr 

loi h 

modified Bessel functions of the first kind 

Pr 

k 

turbulence kinetic energy, m 2 /s 2 , 

Psat 


wave number, 1/m, or 

Q 


thermal conductivity, J/(ms K) 

V 

kij 

binary interaction parameter 

Tk 

Vko 

K b 

a breakup constant, 1/s 

K sb 

wave number associated with 

sheet breakup, 1/m 

rSMR 

R p 

ki,k 2 

constants 

Ry 

K Lb 

wave number associated with 

Re 


ligament breakup, 1/m 

RND 

Dq , K { 

modified Bessel functions of the 

Sh 

l 

second kind 

modified wave number, 1/m 

Sk 

Smlc 

Ik 

mixture latent heat of evaporation, J /kg 

lk,eff 

effective latent heat of evaporation, 

Smle 


J/kg (defined in Eq. (35)) 

Ikn 

heat of vaporization at normal 

Smlm 

L 

boiling point, J/kg 

liquid sheet breakup length, m 


liquid mass flow rate, kg/s 

droplet vaporization rate, kg/s 

droplet vaporization rate under 

flash evaporating conditions, kg/s 

initial mass flow rate associated 

with kth droplet group, kg/s 

droplet vaporization rate due 

to heat transfer, kg/s 

mass of the parent drop, kg 

mean mass of the product drop distribution, kg 

molecular weight, kg/kg-mole 

molecular weight of gas excluding 

fuel species, kg/kg-mole 

molecular weight of fuel, kg/kg-mole 

molecular weight of ith 

species, kg/kg-mole 

shed drop mass, kg 

mass mean diameter, m 

number of droplets in kth group 

drop number in the parent group 

drop number in the product group 

non-dimensional number (=mo/fh(i)) 

drop number 

parent drop number 

number of surfaces contained in 

a given computational cell 

total number of computational cells 

Nusselt number 

pressure, N/m 2 

critical pressure, N/m 2 

Prandtl number 

reduced pressure, P/P c 

saturation pressure, N/m 2 

density ratio {—p g /pi) 

radius (=- v /a 3 /tsmr), m 

droplet radius, m 

initial droplet radial location, m 

Sauter mean radius, m 

switch factor defined in Eq. (61) 

gas constant, J/(kg K) 

Reynolds number 

random number 

Sherwood number 

droplet radius-squared ( = r 2 .), m 2 

source term contribution from liquid 

exchange to mass conservation, kg/s 

source term contribution from liquid 

exchange to energy conservation, J /s 

source term contribution from liquid 

exchange to momentum conservation, kg m/s 2 


vii 
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Smis source term contribution from liquid 
exchange to species conservation, kg/s 

SAID Sauter mean diameter, m 
t time, s 

td non-dimensional time ) 

T temperature, K, or 

non-dimensional number (=ZWe°g ) 

Tf, boiling temperature, K 

T n b normal boiling temperature, K 

T c critical temperature, K 

X/ kth droplet temperature, K 

T r reduced temperature, T /T c 

Ud drop deceleration, m/s 2 

Ui ith velocity component, m/s 

Uik ith velocity component of 

kth droplet group, m/s 

U relative velocity between gas and liquid, 
m/s, or relative drop velocity, m/s 
vt product drop velocity component, m/s 

V initial liquid velocity, m/s 

V c volume of the computational cell, m 3 , or 

critical molar volume, m 3 /kg-mole 
V„ molar volume at normal pressure, 

m 3 /kg- mole 

We Weber number (= pgh J- r ) 

ibj gas-phase chemical reaction rate, 1 /s 

x Cartesian coordinate, m, or drop 

deformation distance, m 

Xi Cartesian coordinate in the ith direction, m 

y Cartesian coordinate, m, or non-dimensional 

deformation parameter (=— ) 
yj mass fraction of jth species 

x spatial vector 

z Cartesian coordinate, m 

Z compressibility factor (=^), or 

non-dimensional number 
Z c critical compressibility factor (= ) 

a represents a coordinate related to 

a Hill’s Vortex, spray cone rotation 
angle, deg., or a constant 

a s overall heat transfer coefficient, kJ/sm 2o K 
ot\ a constant 

«2 a constant 

/ 3 spray cone rotation angle, deg. 

X mole fraction 

S Dirac-delta function, or 

initial liquid sheet thickness, m 

A p pressure drop in the injector, N/m 2 

A tf local time step used in the flow solver, s 


A t g i global time step used in the spray solver, s 

A tu fuel injection time step, s 

At m i allowable time step in the spray solver, s 
AV computational cell volume, m 3 

e rate of turbulence dissipation, m 2 /s 3 

ej fractional mass evaporating rate of species 
at the droplet surface 
77 wave amplitude, m 

turbulent diffusion coefficient, kg/ms, or 

a factor (=210 \ Tc pl \ l ^) 

A thermal conductivity, J / (ms K) , or 
wavelength, m 

A wavelength associated with fi, or 

wavelength associated with 
Rayleigh-Taylor breakup, m 

/i dynamic viscosity, kg/ms 

v kinematic viscosity, m 2 /s 

u turbulence frequency, 1 /s, 

frequency (=^yj - pr), 1 /s, 
complex growth rate, 1 /s, or 
acentric factor 

n maximum growth rate, 1 /s 

p density, kg/m 3 

pin liquid density (at 1 bar, 273.15 K), kg/m 3 
p r reduced density (=p/p c ) 

a potential length constant, Angstrom (=0.1 

surface tension, kg/s 2 , or characteristic 
diameter of a molecule in Table 1 
r stress tensor term, kg/ms 2 , or 

characteristic breakup time (=1 /Kb), s 
6 void fraction, or 

spray cone angle, deg. 

Subscripts 

A A-th species component 
b breakup conditions 

A B-th species component 
c critical conditions 

/ fuel 

g gas or global 

i i-th coordinate, i-th species 

component, summation index, or 
imaginary component 
inj injector 

j j-th coordinate, j-th species 

component, or summation index 
k k-th droplet group, k-th coordinate, 
summation index, or liquid conditions 
associated with k-th droplet group 


viii 
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Superscripts 


l liquid 

L liquid ligament 

in multi-component mixture 

n normal, or nth-face of the 

computational cell 
o initial conditions, 

orifice exit conditions, 
injector initial conditions, or 
oxidizer 

p particle, drop, or conditions 

associated with a grid cell 
r reduced, radial coordinate, or 

real component 
RT Rayleigh-Taylor 

s droplet surface, or adjacent 

computational cell 
t time 

x axial or x-coordinate 

y y-coordinate 

z z-coordinate 

, partial differentiation with respect 

to the variable followed by it 


mean, or average 

first order differentiation, or 

flow rate 

second order differentiation 
// fluctuations 
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LSPRAY-III: A Lagrangian Spray Module 


M.S. Raju 

ASRC Aerospace Corporation 
Glenn Research Center 
Cleveland, Ohio 44135 


1 INTRODUCTION 

There are many occurrences of sprays in a vari- 
ety of industrial and power applications and materi- 
als processing [1], A liquid spray is a two phase flow 
with the gas as the continuous phase and the liquid 
as the dispersed phase in the form of droplets or lig- 
aments [1]. The coupling between the two phases oc- 
curs through the exchanges of mass, momentum, and 
energy involving a wide range of thermal, mass, and 
fluid dynamic factors. A number of finite-difference 
formulations have been advanced over the years for 
predicting the flow (mass and momentum) and ther- 
mal properties of a rapidly vaporizing spray. Some 
of the pros and cons of various formulations can be 
found in [1-5]. 

The state of the art in multi-dimensional com- 
bustor modeling, as evidenced by the level of sophis- 
tication employed in terms of modeling and numeri- 
cal accuracy considerations, is also dictated by the 
available computer memory and turnaround times 
afforded by present-day computers. With the aim 
of advancing the current multi-dimensional com- 
putational tools used in the design of advanced 
technology combustors, a solution procedure is de- 
veloped that combines the novelty of the coupled 
CFD/ spray /sc alar Monte Carlo PDF (Probability 
Density Function) computations on unstructured 
grids with the ability to run on parallel architec- 
tures. In this approach, the mean gas-phase velocity 
and turbulence fields are determined from the solu- 
tion of a conventional CFD method, the scalar fields 
of species and enthalpy from a modeled PDF trans- 
port equation using a Monte Carlo method, and a 
Lagrangian-based dilute spray model is used for the 
liquid-phase representation. 

The gas-turbine combustor flows are often char- 
acterized by a complex interaction between various 
rate-controlling processes associated with turbulent 
transport, mixing, chemical kinetics, evaporation and 
spreading rates of spray, convective and radiative 
heat transfer, among others [10]. The phenomena 
to be modeled as controlled by these processes often 
interact with each other at various disparate time and 
length scales. In particular, turbulence plays an im- 
portant role in determining the rates of mass and heat 
transfer, chemical reactions, and liquid-phase evapo- 
ration in many practical combustion devices. The in- 
fluence of turbulence in a diffusion flame manifests it- 
self in several forms, ranging from the so-called wrin- 
kled or stretched flamelets regime to the distributed 


combustion regime, depending upon how turbulence 
interacts with various flame scales [69-70]. 

Most of the turbulence closure models for reac- 
tive flows have difficulty in treating nonlinear reaction 
rates [69-70]. The use of assumed shape PDF meth- 
ods was found to provide reasonable predictions of 
pattern factors and NO.y emissions at the combustor 
exit [71]. However, their extension to multi-scalar 
chemistry becomes quite intractable. The solution 
procedure based on the modeled joint-composition 
PDF transport equation has an advantage in that it 
treats the nonlinear reaction rates without any ap- 
proximation. This approach holds the promise of 
modeling various important combustion phenomena 
relevant to practical combustion devices such as flame 
extinction and blow-off limits, and unburnt hydrocar- 
bons (UHC), CO, and NOx predictions [71]. 

With the aim of demonstrating the viability of 
the PDF approach to the modeling of practical com- 
bustion flows, we have undertaken the task of extend- 
ing this technique to the modeling of sprays as a part 
of the NCC (National Combustion Code) develop- 
ment program [7-15]. In order to facilitate large-scale 
combustor computations, we have extended our pre- 
vious work on the combined CFD/spray/PDF com- 
putations to parallel computing [10-15]. The use 
of parallel computing offers enormous computational 
power and memory as it can make use of hundreds 
of processors in concert to solve a complex problem. 
The trend towards parallel computing is driven by 
two major developments: the widespread use of dis- 
tributed computing and the recent advancements in 
MPPs (Massively Parallel Processors). The solver is 
designed to be massively parallel and automatically 
scales with the number of available processors. A cur- 
rent status of the use of the parallel computing in tur- 
bulent reacting flows involving sprays, scalar Monte 
Carlo PDF, and unstructured grids was described in 
Ref. [11]. It outlined several numerical techniques 
developed for overcoming some of the high computer 
CPU-time and memory-storage requirements associ- 
ated with the use of Monte Carlo solution methods. 
The parallel performance of both the PDF and CFD 
modules was found to be excellent but the results 
were mixed for the spray computations showing rea- 
sonable performance on massively parallel computers 
like Cray T3D; but its performance was poor on the 
workstation clusters [10]. In order to improve the 
parallel performance of the spray module, two differ- 
ent domain decomposition strategies were developed 
and the results from both strategies were summarized 
[10-14], 
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It is also well known that considerable effort 
usually goes into gridding up of complex gas-turbine 
combustor geometries. In order to allow repre- 
sentation of complex geometries with relative ease, 
we have extended our previous work on the com- 
bined CFD/spray/PDF computations to unstruc- 
tured meshes [11-15]. The grid generation time asso- 
ciated with gridding up practical combustor geome- 
tries, which tend to be very complex in shape and 
configuration, could be reduced considerably by mak- 
ing use of existing unstructured grid generators. The 
solver accommodates the use of an unstructured mesh 
with mixed elements: triangular and/or quadrilateral 
for 2D (two-dimensional) geometries and tetrahedral 
for 3D. 

With the development of two computer mod- 
ules, LSPRAY - a Lagrangian spray solver [14] and 
EUPDF - an Eulerian Monte Carlo PDF solver [15], 
we were able to demonstrate the use of the joint 
scalar Monte Carlo PDF method in the modeling of 
complex multidimensional reacting flows (e.g., gas- 
turbine combustor flows) [10-15]. In this manual, we 
only concentrate on providing complete details of our 
spray module along with some more details on the 
high pressure modifications made to the gas-phase 
calculations. However, further details on the appli- 
cation of the joint scalar Monte Carlo PDF method 
to two-phase flows could be found elsewhere in [10- 
13,15]. 

In this manual, we summarize some important 
aspects of our spray formulation without making any 
attempt to provide an in-depth review on the fluid 
dynamic and transport behavior of reacting sprays 
[6-15]. Depending on the nature of the spray, an 
appropriate selection could be made from the choice 
of various well-known spray formulations (multicon- 
tinua, discrete-particle, or probabilistic) based on ei- 
ther a Lagrangian or an Eulerian representation for 
the liquid-phase equations by making use of appro- 
priate droplet sub-grid models. The present solu- 
tion procedure could be used within the context of 
both multicontinua and probabilistic spray formula- 
tions, as it allows for resolution on a scale greater 
than the average spacing between two neighboring 
droplets [1]. For NCC, the adopted choice for the 
gas phase was an Eulerian scheme. The liquid-phase 
equations form a system of hyperbolic equations and 
they could be solved by means of either an Eulerian or 
a Lagrangian representation. A Lagrangian scheme is 
chosen as it reduces the errors associated with numer- 


ical diffusion. The liquid-phase formulation is based 
on various well-established models for droplet drag; 
the vaporization models of a poly disperse spray take 
into account the transient effects associated with the 
droplet internal heating and the forced convection ef- 
fects associated with droplet internal circulation; and 
it employs models for gas-film valid over a wide range 
of low to intermediate droplet Reynolds numbers [7] . 
And the present stochastic particle tracking method 
is applicable for flows with a dilute spray approxima- 
tion where the droplet loading is low. The numeri- 
cal method could be used within the context of both 
steady and unsteady calculations [8-13]. Not consid- 
ered in the present release of the code are the effects 
associated with droplet/shock interaction and dense 
spray effects. 

Currently, most of the finite-difference tech- 
niques used for predicting the two-phase flows make 
use of the physics derived from single-component liq- 
uid droplet studies with constant properties. How- 
ever, it is well known that most of the gas-turbine fu- 
els are multicomponent mixtures of many compounds 
with a wide distillation curve [7,18]. The multicom- 
ponent nature of the liquid sprays is becoming evi- 
dent with the increasing need to use jet fuels derived 
from heavier petroleum compounds. The gasifica- 
tion behavior of a multicomponent fuel droplet may 
differ significantly over that of a pure single compo- 
nent fuel droplet [18]. Also, the calculation of the 
variable thermo-transport properties of the liquid- 
mixtures becomes more important at high pressures. 
The flame ignition characteristics (such as the phe- 
nomena assciated with flame blow-off and extinction 
conditions) could also be influenced by the nonuni- 
form concentration of the fuels with different volatili- 
ties. However, the importance of the multicomponent 
liquid fuels with variable properties received little at- 
tention in the modeling of comprehensive gas-turbine 
combustor spray computations. With this in mind, 
we have made the following improvements: 

(a) In order to deal with the gas turbine fu- 
els that are mixtures of many compounds, we have 
extended the spray formulation to multicomponent 
liquid sprays. This implementation also takes into 
account the effect of variable liquid properties. 

(b) In order to take into account the nonideal 
gas behavior under high pressure conditions, we have 
completed the integration of the Peng-Robinson EOS 
into our CFD module. Also, we have completed the 
implementation of a high-pressure correction into the 
calculation of transport properties in the gas phase. 
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These modifications would enable the calculation of 
high-pressure flow properties in the gas phase. 

In order to reduce some uncertainty associated 
with the specification of initial droplet conditions, we 
have undertaken the task of integrating an atomiza- 
tion module into our spray solution procedure. The 
atomization module was developed by CFDRC Inc. 
[47] in collaboration with the university of Wisconsin 
(UW) [45-46]. Our computational experience with 
the integration of the atomization module was sum- 
marized in [59-60]. A complete description of all the 
primary atomization as well as the secondary droplet 
breakup models contained in the CFDRC/UW atom- 
ization module is provided in the appendix. 

In some applications involving scramjet and 
ramjet afterburners, liquid fuel may be superheated 
before its injection because the same fuel is often used 
as a coolant [61]. Under some gas-turbine conditions, 
it is estimated that a small fraction of the liquid fuel 
may be released by flash boiling, and there are some 
reported cases of flash related engine performance 
problems in gasoline direct-injection internal combus- 
tion engines [62]. Because of its occurrence in some 
NASA-related critical applications, we have recently 
completed the task of implementing a superheat va- 
porization model into our spray solution procedure. 

Some important aspects of the spray module are 
summarized below: 

• An efficient particle tracking algorithm was de- 
veloped and implemented into the Lagrangian 
spray solver in order to facilitate particle move- 
ment in an unstructured grid of mixed elements. 

• LSPRAY-III is currently coupled with an un- 
structured flow (CFD) solver of NCC [21-23], 
and an Eulerian-based Monte Carlo probability 
density function solver - EUPDF [15], which were 
selected to be as the integral components of the 
NCC cluster of modules. EUPDF rpovides the 
solution for the scalar fields of species and en- 
thalpy from a modeled PDF transport equation 
using a Monte Carlo method. 

• The spray solver receives the mean velocity and 
turbulence fields from the flow solver. The so- 
lution for the scalar (energy and species) fields 
could be provided by means of either a conven- 
tional CFD solver or a Monte Carlo PDF solver 
depending on the choice of the solver. 

• The spray solver supplies the spray source-term 
contributions arising from the exchanges of mass, 


momentum and energy with the liquid-phase 
to the flow solvers (CFD and/or Monte Carlo 
PDF). This information could be used in either 
conservative or non-conservative finite-difference 
formulations of the gas phase equations. 

The furnished code demonstrates the success- 
ful methods used for parallelization and coupling of 
the spray to the flow code. First, complete details of 
the spray solution procedure is presented along with 
several other numerical issues related to the coupling 
between the CFD, LSPRAY-III, and EUPDF solvers. 
It is followed by a brief description of the combined 
parallel performance of the three solvers (CFD, EU- 
PDF, and LSPRAY-III) along with a brief summary 
of the validation cases. 

2 GOVERNING EQUATIONS FOR GAS 
PHASE 

Here, we summarize the governing gas-phase conser- 
vation equations in Eulerian coordinates [1]. This 
is done for the purpose of identifying the interphase 
source terms arising from the exchanges of mass, mo- 
mentum, and energy with the liquid phase. They are 
valid for a dilute spray with a void fraction of the 
gas, 9 , close to unity. The void fraction is defined as 
the ratio of the equivalent volume of gas to a given 
volume of a gas and liquid mixture. 

The conservation of the mass leads to: 

[, pv c ],t + [pVcUi\, Xi = Smic = ^2 n k rh k (1) 

k 

The source term is given as a summation over dif- 
ferent classes of droplets. Each class represents the 
average properties of a different polydisperse group of 
droplets. Here, n k represents the number of droplets 
in a given class and fn k represents the corresponding 
mass vaporization rate. 

For the conservation of the jth species, we have: 

\P V cVj\,t + \pVcUiyj] iXi - [pV c Dy jtXi ] tXi - pV c Wj = 

Smls = ^ ^ Wlk ( 2 ) 

k 

where 

Wj = 0 and Cj = 1 

3 3 
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For the species conservation, the source term con- 
tains an additional variable, Cj, which is defined as 
the fractional vaporization rate for species j. 

For the momentum conservation, we have: 


[pV c ‘Ui],t + [pV c UiUj] tXj + \pV c ] >Xi - 

[$ Xj [(1 b) V c Tnj ] , X j — Smlm 

n k rhk U ki - ^2 y Pk r\ n k u kijt (3) 
k k 

where the shear stress in Eq. (3) is given by: 


T~ij ktj.xi] I ,x :! 

For the momentum conservation, the first source 
term represents the momentum associated with liquid 
fuel vapor and the second represents the momentum 
change associated with droplet drag. 

For the energy conservation, we have: 


[pVch] ,t+[pV cUi h} , Xi - \eV c \T Xi } , Xi - [(1-0) V c \iT Xi \, Xi 

= Smle = ^ ' kl k fll k ^ h s ^k,eff^j (4) 

k 

Similarly, the energy conservation has a source term 
contribution from the added liquid fuel vapor and an 
additional contribution associated with the effective 
latent heat of vaporization which accounts for the 
heat flux (loss or gain) to the droplet interior from 
the ambient. 

The main purpose of the spray solver is to cal- 
culate the source terms arising from the exchanges 
of the mass, momentum, and energy. In the case of 
NCC, it provides the calculated source terms to both 
the CFD and Monte Carlo PDF solvers. 

3 HIGH PRESSURE EQUATION OF 
STATE 

In order to calculate the high pressure gas be- 
havior, the Peng-Robinson EOS (Equation-Of-State) 
is employed for a multi-component mixture in the fol- 
lowing form [16-17,28]: 


P = 


RT 


V-b 


m 


V 2 + 2 b m V - b 2 m 


where 


( 5 ) 


a m = E, Ej yi‘yj( a i a j) 1/2 ( l - hj), 

bm — E? y* ' 




Qi = 


0.07780 RT ic 

Pic ’ 

0.45724i? 2 T? 


[1 + fU 1 - El 72 )] 2 , 


Tir = T /T ic , 


f VjJ = 0.37464 + 1.54226^ - 0.26992 Wi 2 , 

u>i is known as the acentric factor of the molecules 
which is a measure of non-sphericity of the molecules, 
and kij is known as the binary interaction coefficient. 
However, the Peng-Robinson EOS is rewritten as a 
cubic EOS in terms of the compressibility factor, Z 
(= ^), before it is solved: 

Z z - (1 - B*)Z 2 + (A* - 2B* - 3B* 2 )Z 

- A*B * + B* 2 + B* 3 = 0 (6) 

where 

A* = pT, and B* = b -%f. 

We have chosen the Peng-Robinson EOS be- 
cause of its simplicity and, more importantly, it 
proved to be very useful in the supercritical droplet 
vaporization studies of [19-20]. 

Table 1 lists various physical constants for some 
of the species that we found to be of interest in our 
spray computations. It contains values for the boiling 
temperature at normal pressure, the critical tempera- 
ture, the critical pressure, the critical density, the la- 
tent heat of vaporization at normal pressure, the crit- 
ical molar volume, the molar volume at normal pres- 
sure, and the acentric factor of the molecules. Most 
of this data is collected from [16-18] and is useful in 
both the evaluation of the Peng-Robinson EOS and 
the calculation of various other variable properties. 
And Table 2 lists the binary interaction parameters, 
kij, used in a calculation for the multi-component 
mixture of n-heptane, 0 2 , N 2 , C0 2 , & H 2 0. While 
the data for most of the binary pairs was obtained 
from various reference books, some of the missing 
data, however, is replaced with known data found for 
other binary pairs of molecules with similar molecular 
weights. 
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Table 1. Physical constants. 

Species 

T n b 

T c 

Pc 

Pc 

^kn 

K c 

v n 

LO 

a 

a / k 


(K) 

(K) 

(atm) 

(Kg/m 3 ) 

(KJ/kg) 

(cm 3 /g-mole) 

(cm 3 /g-mole) 


(A °) 

(K) 

CqHi4 

341.9 

507.4 

30.0 

660.0 

334.8 

370.0 

140.06 

0.296 

5.949 

399.3 

C-jHiq 

371.6 

540.2 

27.0 

682.0 

316.3 

432.0 

162.00 

0.351 

6.297 

419.031 

CsHis 

398.82 

568.8 

24.6 

718.5 

301.3 

492.0 

188.8 

0.394 

6.62 

488.15 

C10H22 

447.3 

617.6 

20.8 

728.3 

276.1 

603.0 

233.68 

0.490 

7.16 

540.06 

Cl 2^26 

489.5 

658.3 

18.0 

748.0 

256.3 

713.0 

278.54 

0.562 

7.655 

583.68 

C14H30 

526.7 

694.0 

16.0 

763.0 

240.1 

830.0 

326.62 

0.679 

8.067 

629.08 

n 2 

77.4 

126.2 

33.9 

807.1 

197.6 

90.1 

31.87 

0.039 

3.681 

91.5 

o 2 

90.2 

154.6 

50.4 

1135.7 

212.3 

74.4 

26.08 

0.025 

3.433 

113.0 

C 02 

00.0 

304.1 

73.8 

000.0 

000.0 

94.0 

33.32 

0.239 

3.996 

190.0 

h 2 o 

373.2 

647.3 

221.2 

958.1 

2257.2 

57.1 

19.76 

0.344 

2.641 

809.1 


Table 2. Binary Interaction Parameters, 


C7H1Q 

0=1) 

0 2 

0=2) 

n 2 

0=3) 

C0 2 

0=4) 

h 2 o 

0=5) 

C 7 H \ 6 

(i=l) 

0.0000 

0.1321 

0.1440 

0.1000 

0.1484 

O 2 

(i=2) 

0.1321 

0.0000 

-0.0119 

-0.0289 

0.0910 

n 2 

(i=3) 

0.1440 

-0.0119 

0.0000 

-0.0170 

0.1030 

C0 2 

(i=4) 

0.1000 

-0.0289 

-0.0170 

0.0000 

0.1200 

h 2 o 

(i=5) 

0.1484 

0.0910 

0.1030 

0.1200 

0.0000 
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4 HIGH PRESSURE CORRECTIONS FOR 
GAS-PHASE TRANSPORT PROPERTIES 


The effect of pressure on the viscosity of pure 
gases is determined by making use of the Reichenberg 
method [16,19,29]: 


P^ 

/^n 


1 + Q v 


A„P, 


, 3/2 


B v P r + (1 + C„P r °”) 


*M-i 


( 7 ) 


where 


A v = ^exp(a 2 T“) B v = A v {^T r - fa) 
C v = %-expinT?) D v = ^exp(A 2 T r d ) 


ai = 1.982510 -03 
/3i = 1.6552 
7i = 0.1319 
Ai = 2.9496 


a 2 = 5.2683 
(3 2 = 1.2760 
y 2 = 3.7035 
A 2 = 2.9190 


a = -0.5767 

c = -79.8678 
d= -16.6169 



0.0 1.0 2.0 3.0 4.0 5.0 6.0 

Reduced pressure, P r 


Figure 1 Takahashi correlation showing the 
variation of the binary diffusion coefficient 
versus reduced pressure at different reduced 
temperatures. 


and Q v = 1.0 for non-polar molecules. T r 

The effect of pressure on the thermal conduc- T c ab 
tivity of pure gases is determined by making use of P r 
the Stiel & Thodos method [16,19,30]: P c ab 


T 


T cA b 

UaT c a + UbT c b 

p 


P cAB 

VaPcA + UbPcB 


(A - \ n )TZl = 1.2210 -2 [exp(0.535/9 r ) - 1] (8) 

when p r < 0.5 

(A - A„)rZ c 5 = 1.1410 -2 [exp(0.67p r ) - 1.069] (9) 
when 0.5 < p r < 2.0 


(A - \ n )TZl = 2.6010 _2 [exp(1.155p r ) + 2.016] (10) 
when 2.0 < p r < 2.8 


where A is in W/(m.K), T = 210 [ T p/ ] 1 ^ 6 , Z c is the 
critical compressibility, and p r is the reduced density 

p/Pc = V c /v. 

The effect of pressure & temperature on diffu- 
sion coefficients in gases is determined by making use 
of Takahashi correlation [16,19,31]: 


DabP 

(■ D ab P)+ 


f(Tr,Pr) 


( 11 ) 


where Dab = diffusion coefficient, cm 2 /s, P = 
pressure, bar, the superscript + indicates that low- 
pressure values are to be used, and the reduced pres- 
sure and temperature in the above equation are cal- 
culated as follows: 


This correlation is valid for a system involving 
a trace solute diffusing in a supercritical fluid. It is 
shown to provide satisfactory results but it is based 
on a very limited set of available experimental data. 
The graphical representation of the Takahashi func- 
tion is given in Fig. 1. 

It is noteworthy that at low pressures, the poly- 
nomial fits for the variable thermodynamic properties 
are taken from the data set compiled by McBride et el 
[24]. The transport properties involving the thermal 
conductivity and molecular viscosity for individual 
species is estimated based on the Chapman-Enskog 
collision theory [25-27]. And Wilke’s formulae is used 
to determine the properties of mixture [25-27]. The 
binary-diffusion coefficients are determined based on 
the Chapman-Enskog theory and the Lennard-Jones 
potential [25-27]. 

5 LIQUID-PHASE EQUATIONS 

Here, we summarize the governing equations for 
the liquid-phase based on a Lagrangian formulation 
where the equations for particle position and velocity 
are described by a set of ordinary differential equa- 
tions. For the particle position of the ktlr droplet 
group, we have: 
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dx ik 

dt. 

For the droplet velocity: 


— n ik 


duik 3 G jg pg S Re k 
dt 


16 


Pkr 2 k 


Wig Ug V. ; k 


where 


24 


,_.. h Ref' 
Cd ^r^{ 1 + — 


(12) 


(13) 


(14) 


Rek = 2 


r_kPg_ 

Pgs 


[(u g + Ug - Uk ) . ( Ug + Ug - Uk)] 


1/2 


Pk ^ ) VkiPki 


(18) 


and the individual component liquid density is given 
by: 


Pki — 


Mki 

vi~ 


where the molar volume, Vki, is 


V ki = V ci { 0.29056 - 0.8775 w 4 ) c 


(19) 


(20) 


According to Yuen and Chen [32] , the droplet drag for 
a vaporizing droplet could be calculated by a solid- 
sphere drag correlation but suggested using a correc- 
tion in the evaluation of the droplet Reynolds number 
for the average viscosity based on a 1/3 weighting rule 
as given by Eq. (46). The droplet Reynolds number 
is given by: 


and 


C v 


T k 

1 - 

T ■ J 

r.i. J 


(15) 

By following the approach taken from KIVA 
[33], a gas turbulence velocity, u' g , is added to the 
local mean gas velocity when calculating the droplet 
drag and vaporization rates. The gas velocity fluctu- 
ations is calculated by randomly sampling a Gaussian 
distribution with mean square deviation, 2/3 k. The 
Gaussian is given by 


G(u' g ) = (4/3nk) 3 / 2 exp[— 3\u' g \ 2 /4k\ (16) 

The gas fluctuating component is calculated once at 
every turbulence interaction time, f tur , and is oth- 
erwise held constant [33]. The correlation time is 
taken to be the minimum of either the eddy time or 
the transit time taken by the droplet to traverse the 
eddy. It is given by 


The droplet regression rate is determined from 
one of three different correlations depending upon the 
droplet Reynolds number range. The first correlation 
as given by Eq. (21) is based on a gas-film analysis 
developed by Tong & Sirignano [34]. It is based on a 
a combination of stagnation and flat-plate boundary- 
layer analysis and is valid for Reynolds numbers in 
the intermediate range. The last two correlations as 
given by Eqs. (22) & (23) are taken from Clift et 
al [35] and are valid when Rek < 20. In fact, when 
the droplet Reynolds number goes to zero, Eq. (23) 
becomes identical to the droplet regression rate of a 
vaporizing droplet in a quiescent medium [27] . 


ds k 

dt 


= -2 


PgBg 

Pk 


1 1/2 


Re b 


f{B k ) (21) 


if Rek > 20 


dsk 

dt 


PgDg 

Pk 


1 + (1 + - Refc) 1 / 3 


Rel 077 


ln( 1 


if 1 < Rek < 20 


B k ) 

(22) 


f h 3 / 2 1 

t-tur min ( , Ctt -j ; “ ; 

Ve e \Ug + u'-Uk\ 


(17) 


where Ctt. is an empirical constant with a value of 
0.16432. 

The liquid mixture density, p k , in Eq. (13) is 
given by [16,18]: 


dsk 

dt 


PgBg 


Pk 

if Re k < 1 


1 + (1 + i?e fc ) 1/3 l ln{l + B k ) (23) 


where B k is the Spalding mass transfer number and 
is given by: 
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Bk 


( Vfs - Vf) 
(! - Vfs) 


(24) 


and yf s is given as a summation over the fuel-species 
mass fractions at the droplet interface: 


yfs = y « ( 25 ) 

i 

and yf is given as a summation over the fuel-species 
mass fractions in the ambient: 


Vf = Yl y B ( 26 ) 

i 

and the function f(Bk) is similar to that of a Blassius 
function [1, 36] and is obtained from an analysis simi- 
lar to that of Emmon’s boundary-layer flow over a flat 
plate with blowing [36]. The range of validity of this 
function was extended in Raju and Sirignano [7] to 
consider the effects associated with a boundary-layer 
flow with suction. 

The internal droplet temperature is determined 
based on a vortex model [34] . The governing equation 
for the internal droplet temperature is given by: 


dT k A; 

dt Cpipirl 


^ dTk 

a „ 0 + (1 + C(t)a) 


da 2 


da 


where 


CM = ^ 


C pW i 

A / 


(27) 


dTk 

(28) 


where a represents the coordinate normal to the 
stream-surface of a Hill’s Vortex in the circulating 
fluid, and C{t) represents a nondimension al form of 
the droplet regression rate. The initial and boundary 
conditions for Eq. (27) are given by: 


t — t injection i Bk — 7 k . c 


a = 0, 


dT k 

da 


1 

17 


Cpipi 

Ai 


,dT k 
: dt 


(29) 

(30) 


dT k _ 3 pk r . , i ds k 

h da ~ 32 A, [h ’ eff h] dt (31) 


and the individual component latent heat of vapor- 
ization, Iki, is given by [16,18]: 


l 


ki 


/ T c j — Tfc 

\ r r 'T' 

v -*■ ci - L bi 


0.38 


(33) 


and the droplet boiling temperature is given by 


rp IkinMj/ Ru 

M ~ h,„M, - HP) 1 j 

and, finally, the effective latent heat of vaporization, 
lk,effi is defined as: 

h,eff=lk+ (35) 

m k V dr J s 

It is a very useful parameter as it represents the total 
energy loss associated with the latent heat of vapor- 
ization in addition to the the heat loss to the droplet 
interior. h,eff is calculated by means of the following 
relationship [18]: 


, C p (T g ~T ks ) 

k ’ eff (1 + B k y/ Le - 1 


(36) 


Similar to the internal thermal transport, the 
internal mass transport of a multi-component fuel is 
given by: 


dyki 

dt 



a 


d 2 y ki 
da 2 


+ (1 + C(t)a)^ 


(37) 


The initial and boundary conditions for Eq. (37) are 
given by: 


t — t i n j ec tiom II k i y k i,o 


a = 0, 


dy k i 

da 


1 

17 


D k 


dyki 

dt 


(38) 

(39) 


1 d d ki 3 1 r l d3k OM 

° = 1 ' = (40) 


By knowing the mass fractions of the liquid 
species at the droplet surface, the corresponding mole 
fractions are determined by 


where a = 0 refers to the vortex center, and a = 1 
refers to the droplet surface, and the mixture latent 
heat of vaporization, l k , is given by 

h = ^ f-ihi (32) 

i 


%iks 


Viks / ’M-i 
X)?: Vyks/Mi 


(41) 


At the droplet interface, the mole fractions of 
the gas species are obtained by means of Raoult’s 
law: 
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%is — p^iks^is (42) 

where the partial pressure, P,; s , is determined by 
means of the Clausius-Clapeyron relationship: 


P is = exp 


Iki j 

( i 

_ j_Y 

_Ru 

\T bl 

Tks J _ 


(43) 


Then, the corresponding mass fractions in the gas- 
phase at the droplet interface are given by 


v . = XiaMi (44) 

M a ( 1 - x is) + M iXis 

where M a is the molecular weight of the gas excluding 
fuel vapor. And the fractional mass vaporization rate 
of liquid species, e,, is given by 

£i = Vis + (1 ^ Ufa) — — — (45) 

Vfs - VS 

It is noteworthy that the thermodynamic and 
transport properties at the gas film are calculated at 
the temperature and composition as determined by 
the following one-third rule: 

1 2 

fyavg = 4“ (46) 

The correlations for the gas-phase thermody- 
namic and transport properties are described in Sec- 
tions 2 & 4. In a similar way, the liquid-phase ther- 
modynamic and transport properties are determined 
based on the correlations described in the next sec- 
tion. 


6 LOW PRESSURE LIQUID 
THERMODYNAMIC & TRANSPORT 
PROPERTIES 

The specific heat at constant pressure, C p i , 
thermal conductivity, A i, and viscosity, /.q, are evalu- 
ated by means of the following expressions: 


C p i 

II 

£ 

O 

+ CpiiT + C p i2 

T 2 + 

CpisT 3 + 

C p uT a 







(47) 

A i = 

= A/o + 

A;i T + X l2 T 2 + 

a, 3 t 3 

+ 

A m T 4 + A l5 T 5 







(48) 


In m = 

= Mzo + Hn/T + 

P 12 T 

+ 

mT 2 

(49) 


where C p i is in J/(kg K), pi in (/xPA s), and A; in 
(W/m K). 

Tables 3-5 provide the polynomial constants 
used in Eqs. (47)-(49) for some of the species that 
we found to be of interest in our spray computations. 
Table 3 provides the constants for the liquid specific 
heat, C p i, Table 4 for the liquid thermal conductivity, 
A i, and Table 5 for the liquid molecular viscosity, /./; . 
These tables are compiled with the data taken mostly 
from the references of [16,18]. 

The binary diffusion coefficient, Dij, is evalu- 
ated as follows [16]: 


— 


K dif T 


(50) 


where K di j = 8.210 _s (l + [qy 2 -] 2 ^ 3 )- One should be 
careful in using this approximation as it is based on 
a scarce set of available experimental data. 

The specific heat for a multicomponent mixture 
is given by: 


Cpm = ViCpi (51) 

i= 1 

And the thermal conductivity of a multicomponent 
mixture is calculated by means of the Li method [16]. 

n n 

A m ~ ^ ( y ( (52) 

i= 1 9 = 1 

where 

A ij = 2(A i 1 + A- 1 ) _1 

A. - *iVi 


where x t is the mole fraction of the species i, (j)i is a 
volume fraction of the ith species, and Vj is the molar 
volume of the pure fluid. 

7 SUPERHEAT LIQUID FUEL 
VAPORIZATION 

Flashing phenomena refers to a process that is 
in thermodynamic non-equilibrium when a liquid is 
superheated [72-73]. The reasons for its occurrence 
are mainly two-fold [72-73]: (1) a liquid fuel can be 
heated to a higher temperature (above its saturation 
temperature) while its pressure is maintained, and 
(2) when a liquid is depressurized rapidly it can lead 
to flash injection because the thermal inertia initially 
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Table 3. Polynomial constants for liquid specific heat. 

Fuel 

CpZO 

CpZl 

Cpl2 

Cpl 3 

CpZ 4 

C(>H\4 

2.4169 

-5.9866e-03 

2.0959e-05 

-8.4489e-09 

0.0 

C7H1Q 

4.8227 

-3.6980e-02 

1.6777e-04 

-3.0987e-07 

2.2081e-10 

CgHig 

9.2189 

-8.8314e-02 

3.8869e-04 

-7.2539e-07 

5.0776e-10 

C 10 H 22 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

C 12 H 26 

4.7900 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

Cl4^30 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 



Table 4 

. Polynomial constants for liquid thermal conductivity. 


Fuel 

Azo 

Azi 

Az2 

A/3 

A/4 

A/5 

CqH\4 

0.37078 

-5.4313e-03 

4.628e-05 

-1.8002e-07 

3.2243e-10 

-2.1832e-13 

C7H1Q 

0.13236 

9.4441e-04 

-6.588d-06 

1.4617e-08 

-1.1244e-ll 

0.0 

CsHis 

0.25652 

-7.5401e-04 

1.5872e-06 

-1.6795e-09 

-1.3375e-16 

0.0 

C 10 H 22 

0.22179 

-2.3699e-04 

-6.94e-07 

2.0415e-09 

-1.5741e-12 

0.0 

C 12 H 26 

0.17609 

4.2463e-05 

-7.4467e-07 

6.9446e-10 

0.0 

0.0 

C 14 H 30 

0.18801 

-9.1399e-05 

-2.1464e-07 

1.1655e-10 

0.0 

0.0 

n 2 

-2.629E-1 

-1.545E-3 

-9.450E-7 

0.0 

0.0 

0.0 

0 2 

2.444E-1 

-8.813E-4 

-2.023E-6 

0.0 

0.0 

0.0 

C0 2 

4.070E-1 

-8.438E-4 

-9.626E-7 

0.0 

0.0 

0.0 

h 2 o 

-3.838E-1 

5.254E-3 

-6.369E-6 

0.0 

0.0 

0.0 


Table 5. 

Polynomial constants for liquid molecular viscosity. 

Fuel 



H 12 


CqH \4 

-4.034E+00 

8.354E+02 

0.0000000 

0.0000000 

C 7 H 1 Q 

-4.325E+00 

1.006E+03 

0.0000000 

0.0000000 

CgHig 

-4.333E+00 

1.091E+03 

0.0000000 

0.0000000 

C 10 H 22 

-4.460E+00 

1.286E+03 

0.0000000 

0.0000000 

C 12 H 26 

-4.562E+00 

1.454E+03 

0.0000000 

0.0000000 

C 14 H 30 

-4.615E+00 

1.588E+03 

0.0000000 

0.0000000 

n 2 

-2.795E+01 

8.660E+02 

2.763E-01 

-1.084E-03 

0 2 

-4.771E+00 

2.146E+02 

1.389E+02 

-6.255E-05 

C0 2 

-3.097E+00 

4.886E+01 

2.381E-02 

-7.840E-05 

h 2 o 

-2.471E+01 

4.209E+03 

4.527E-02 

-3.376E-05 
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tends to maintain its bulk internal temperature above 
the saturation temperature. Although flash evapora- 
tion is considered to be detrimental to engine perfor- 
mance under normal circumstances, it can have some 
potential benefits: it is known to produce a fine spray 
with enhanced atomization, increase effective spray 
cone angle, and decrease spray penetration [61]. 

An understanding of flash injection is of im- 
portance in some applications involving scramjet and 
ramjet afterburners because the same liquid fuel is of- 
ten used as a coolant coupled with engine conditions 
where nozzles operate at low back pressures and su- 
personic outflow [61]. The objective of our work is 
to establish a baseline accuracy for existing atomiza- 
tion and vaporization models valid under superheat 
conditions by undertaking a critical review of exist- 
ing experimental data for validation. This work is 
funded through the NASA’s fundamental aeronau- 
tics/ supersonic initiative: high altitude emissions. 

We have started out initially with a modeling 
approach based on the existing models developed for 
superheat vaporization. It is based on an extension 
of the classical Zl 2 -theory, and was adopted from the 
papers of Zuo, Gomes, and Rutland [62] and Schmehl 
and Steelant [63-64]. In the classical evaporation 
model, the thermal energy needed for evaporation is 
mostly furnished by the external heat transfer from 
the surrounding gas. Under superheat conditions, the 
characteristic vaporization time resulting from the 
external heat transfer from the surrounding gas is of 
the same order of magnitude as that resulting from 
the flash evaporation. The energy needed for vapor- 
ization at the droplet surface is partly provided by 
the superheat energy stored within the droplet and 
it is controlled by the droplet internal heat trans- 
fer. The modeling approach differs from the classical 
droplet vaporization models in three important ways: 
(1) under superheat conditions, the droplet surface 
mass fraction, Yf s , approaches unity as the droplet 
surface temperature is maintained at the correspond- 
ing liquid boiling temperature; (2) under superheat 
conditions, all the external heat transfer from the sur- 
rounding gas is made available to the vaporization 
process with no apparent increase in the droplet sur- 
face temperature; and (3) the flow of fuel vapor im- 
parted by flash vaporization partly counterbalances 
the flow generated by external heat transfer and may 
significantly reduce the energy transferred from the 
surrounding gas. 

7.1 Vaporization Model Valid Under Superheat 


Conditions 

Based on the governing equations of conserva- 
tion for an isolated spherically symmetric droplet, 
Zuo et al [62] and Schmehl and Steelant [63-64] 
showed that the total evaporation rate, fh k , can be 
calculated as 


rn k = m k , flash + m k ,t (53) 

where the flash boiled vaporization rate, m k jiash, is 
given by 


rnk, flash = 47r rla 8 ^ k — (54) 

Ik 

where T k is the internal droplet temperature and the 
overall heat transfer coefficient, a s (= kJ/s m 2 °K) 
is given by the Adachi correlation [65]: 

= 0.76(7), - T b ) 0 2e (0 < T ks -T b < 5) 

a a = 0.027(7). - 7)) 2 ' 33 (5 < T ks - T b < 25) (55) 

= 13.8(7) - Ti) 0 ' 39 ( T ks -T b > 25) 


The Adachi correlation is valid over a wide range of 
superheat conditions. 

The vaporization rate due to external heat 
transfer, m k ,t, in Eq. (53) is given by 


k Nu r „ mt flash ■. , 

m k ,t = 2irr fc — -r ini + (1 4 ) )B t ] (56) 

Cp 1+ m k, flash m kt 


where the Spalding heat transfer number, B t , and the 
Nusselt number, Nu, are given by 


Bt = 


C p {T g ^T ks ) 


(57) 


h,eff 

Nu = 2(1 + 0.3 Re 1/2 Prl /3 ) (58) 

The corresponding droplet regression rate, 
is given by 


ds k _ m k 
dt 2 tt r k pi 


(59) 


This model is valid over an entire range of su- 
perheat conditions as long as there is some amount of 
superheat energy available within the droplet (T k > 
T b ). 


7.2 Combined Superheat-Classical Vaporization 
Model 
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Under moderate initial superheat conditions, 
only a fraction of the vaporization takes place under 
superheat conditions ( T k > T b ) and the remainder 
takes place under more stable (non-superheat) evap- 
orating conditions (T k < T b ). So there is a need to 
revert back to a vaporization model valid under sta- 
ble evaporating conditions when the internal droplet 
temperature approaches the boiling temperature. In 
the present calculations, the vaporization rate under 
normal evaporating conditions is evaluated by means 
of a simplified classical Z? 2 -theory: 


m fc = 2nr k p g Df gs Sh ln( 1 + B k ) (60) 

where the Spalding mass transfer number, B k , is 
given by Eq. (24) and the Sherwood number, Sh, 
is given by 

Sh = 2(1 + 0.3 Re 1/2 Sc 1 g /3 ) (61) 

7.3 Internal Droplet Temperature Calculation 

Our experience with the validation studies 
showed us that there is a definite need to include 
a calculation involving the internal droplet tempera- 
ture valid under both superheat and normal evapo- 
rating conditions. In our present calculations, it was 
evaluated by means of a simple infinite conductivity 
model. 


dT k 

dt 


dT k 

dt 


3[^fc,e// Ik] ds k 

2C p irl dt 
if T k < T b , and 


3a s 

f'l.pi (dpi 


(T k - T b ) 


if Tk > T b 


(62) 


(63) 


8 DETAILS OF ATOMIZATION 
MODELING & LIQUID JET/SHEET 
BREAKUP 

The success of any spray model depends a great 
deal on the specification of the appropriate injector 
exit conditions. So far in our spray computations, 
we have relied upon either known experimental data 
or data generated from widely-used correlations in 
specifying the initial droplet conditions. In order to 
reduce the uncertainty associated with the specifica- 
tion of the initial droplet conditions, we have under- 
taken the task of integrating an atomization module 


into our spray calculation procedure. The atomiza- 
tion module is based on some recent progress made 
into the modeling of the atomization process [41-48] . 

The atomization process can be broadly clas- 
sified into two breakup regimes: (1) Within the in- 
jector, the inner-nozzle disturbances due to cavita- 
tion may lead to the formation of a fragmented liq- 
uid, and (2) Once the liquid jet emerges outside of 
the nozzle, it becomes unstable under the influence 
of hydrodynamic instabilities created by the exter- 
nal aerodynamic liquid-gas interaction and it then 
breakups up into ligaments, fragments and droplets 
[43-44]. The widely known hydrodynamic instabil- 
ities are of Rayleigh-Taylor and Kelvin-Helmholtz 
kind [41, 43-44]. The Rayleigh-Taylor instability is 
due to inertia of the denser fluid opposing the system 
acceleration in a direction perpendicular to the in- 
terface of the denser fluid and the Kelvin-Helmholtz 
instability is caused by the viscous forces due to the 
relative motion of the fluids [43-44]. The liquid jet 
breakup is also influenced by the nozzle geometry and 
the thermo-physical properties of the fluid. The ap- 
proach taken here makes use of an extensive analysis 
developed over the years in understanding various as- 
pects of the external aerodynamic liquid-gas interac- 
tion but relies primarily on some known correlations 
when it comes to describing the internal- injector pro- 
cesses [41]. 

Based on a linear instability analysis of a 2D 
viscous incompressible fluid moving thorough an in- 
viscid incompressible gas, Reitz and Bracco [41] char- 
acterized the break-up regimes to be four-fold: (1) 
Rayleigh breakup, (2) first wind-induced breakup, 
(3) second wind-induced breakup, (4) atomization. 
In the first two regimes, drops of sizes greater than 
or equal to the nozzle diameter are produced at dis- 
tances far from the nozzle exit. The last two regimes 
are more important to the atomization studies of our 
interest and drops of sizes smaller than the nozzle di- 
ameter are produced near the nozzle exit. The knowl- 
edge gained from the instability analysis of various 
kinds [41, 45-46] is combined with some experimen- 
tal observations to form the basis for the atomization 
and droplet breakup modeling of liquid sprays. As 
various processes associated with atomization such 
as the jet breakup and the breakup of drops remain 
relatively indistinguishable within a dense spray, the 
jet breakup is modeled by making use of a drop repre- 
sentation approach in which discrete parcels of liquid 
are injected in the form of blobs with a characteris- 
tic size representative of the nozzle diameter instead 
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of tracking an actual intact liquid core at the nozzle 
exit. In the case of a planar or conical liquid sheet, 
the discrete parcels essentially represent liquid liga- 
ments. Before the jet breakup the discrete parcels 
stay inside of the liquid core or sheet but after the jet 
breakup they move independently. The breakup cri- 
terion, atomization rate, drop size and velocity and 
the location of the newly formed droplets are pri- 
marily determined based on an instability analysis 
derived from the conservation equations of mass, mo- 
mentum and energy. The analysis of the jet or sheet 
breakup into ligaments or droplets, the stripping of 
the liquid into fragments or droplets, and the for- 
mation of smaller droplets from further breakup of 
ligaments or fragments are all described under the 
classification of primary atomization. 

Some of the large droplets that are formed im- 
mediately after the primary liquid jet breakup may 
further breakup into smaller droplets under the influ- 
ence of aerodynamic instabilities. The large droplets 
first tend to flatten under the influence of aerody- 
namic pressure. Then large amplitude long wave- 
length waves caused by drop deceleration induce 
a Rayleigh-Taylor instability on the flattened drop 
causing it further to breakup into several relatively 
large-size product droplets. While at the same time 
short surface waves induce a Kelvin-Helmholtz insta- 
bility on the windward side of the parent drop re- 
sulting in the generation of much smaller product 
droplets. The breakup of the larger droplets into 
smaller droplets is described under the classification 
of secondary droplet breakup. 

The CFDRC atomization module contains the 
following four primary atomization models: (1) the 
sheet breakup primary atomization model, (2) the 
blob jet primary atomization model, (3) the air- 
blast atomization model & (4) the BLS (Boundary- 
Layer Stripping) primary atomization model; and the 
following three secondary droplet breakup models: 
(1) the Rayleigh-Taylor secondary droplet breakup 
model, (2) the TAB (Taylor Analogy Breakup) sec- 
ondary droplet breakup model, and (3) the ETAB 
(Enhanced TAB) secondary droplet breakup model. 
The choice between various models depends on the 
specific application. More details of all the models 
contained in the atomization module are described in 
the appendix. 

8.1 The Effect of Flash Vaporization on Primary 
Atomization 

Flash evaporation is known to produce a fine 


spray with enhanced atomization, increase effective 
spray cone angle, and decrease spray penetration [61]. 
Here, we consider the effects of flash evaporation on 
the sheet breakup primary atomization model by fol- 
lowing the approach of Zuo et al [62]. Its effect on 
the initial droplet size generated immediately after 
the first ligament breakup, di 3 , is given as a function 
of both engine pressure and a superheat parameter 
as follows: 

d,s = *„(^)°- 27 [i-x(^)°-“ B ] 0<X(^)°' 135 <1 

atm d f 

( 64 ) 

where di n is the corresponding droplet droplet size 
under normal evaporating conditions without flash 
evaporation, and y is defined as a superheat param- 
eter as follows: 


I{T k ) - I{T b ) 

TO 


(65) 


where / is the internal energy of the liquid. Its value 
varies between 0 < % < 1 with y = 0 referring to zero 
flash evaporation and y = 1 to full flash evaporation. 

In Eq. (64), the increase in di S due to an in- 
crease in engine pressure by a factor of ( p p ) 0,27 
is based on an experimental correlation obtained 
from Lefebvre [66]. It reflects the influence of cham- 
ber pressure on wave propagation as it damps wave 
growth. But the decrease by a factor of (1 — 
y( Pa p m ) 0135 ) is due to a significant reduction in 
droplet size caused by both cavitation and bubble 
growth under flash evaporation conditions. It was 
introduced based on the experimental data obtained 
from VanDerWege et al [67] and Reitz [68]. 

As the liquid approaches boiling, it also causes a 
substantial decrease in both intact liquid core length 
and core droplet size leading to a modification of the 
nominal cone angle, 9 , as given by 


9 = 9 n + (144-0„)y 2 (66) 

where 9 n is in degrees for a spray vaporizing under 
normal conditions without flash evaporation. This 
correlation was developed by based on the experimen- 
tal data of Reitz [68]. Both these correlations were 
developed in conjunction with the sheet breakup pri- 
mary atomization model but their validity with other 
primary atomization models needs further investiga- 
tion. 


NASA/CR— 2008-215290 


13 



9 DESCRIPTION OF INITIAL SPRAY 
CONDITIONS 


The spray computations facilitate the use of 
multiple fuel injectors. The same or a different type 
of liquid fuel can be specified for each one of different 
injectors. The initial droplet temperature is assumed 
to be the same for all different droplet groups of any 
given injector. The liquid fuel injection is simulated 
by introducing a number of discretized parcels of liq- 
uid mass at the beginning of every fuel-injection time 
step, A tu. The following three variables play an im- 
portant role in simulating the injector initial condi- 
tions: 

• no_of_holes() - the number of holes per injector, 

• no_of_streams() - the number of droplet streams 
per hole - it is introduced to distinguish the 
initial velocity variation within different droplet 
classes arising from the geometric considerations 
of a chosen spray, & 

• no_of_droplet_groups() - the number of droplet 
groups per stream - it is introduced to distin- 
guish the droplet-size variation within different 
droplet classes of a polydisperse spray. 


three available options: (1) by providing a complete 
specification of the initial conditions through a spray 
table, (2) by invoking certain pre-defined correlations, 
or (3) by chosing one of four available primary at- 
omization models. When Option (3) is selected, the 
logical variable, atomization(), is set to .true, or oth- 
erwise it is set to .false.. 

Option (1): 

When a spray table is defined, one should pro- 
vide the information on the (x,y z) components of 
the initial droplet location, the (u, v, w) compo- 
nents of the droplet velocity, and the initial mass flow 
rate associated with each different droplet group as 
described in ncc _spray .table. 1 of Table 7. This ta- 
ble provides a complete description of the following 
variables: nos(nJ), (ni,xx jnj, yyjnj, zz_inj, uu_inj, 
vvdnj, wwdnj, rdnj, flcbd), ,ni=l,nos(n_i)). The ini- 
tial inputs specified through a spray table should be 
representative of the integrated averages of the ex- 
perimental conditions [10-13]. 

Option (2): 


The total number of droplet groups introduced at 
the beginning of every different injection time step 
for a given injector is therefore equal to a value 
given by the multiplication of these three vari- 
ables. Further details on some of the spray in- 
put variables can be found in the spray input file, 
ncc jnjector.in.l of Table 6, where the integer num- 
ber followed by the last dot represents the in- 
jector number, which in this case happens to be 
one. The table provides a complete description of 
the following input variables: nolc(n_i), out_string, 
(ymki(n_i,n_l), n_l=l,nolc(n_i)), tdrop(nj), atomiza- 
tion(nJ), drop -breakup _model(nJ), spray _table(n _i), 
steady _spray -model, no_of_holes(nj), 

no_of_streams (n_i) , 

no_of_droplet -groups (n j), lmdis(nj), smdm(ni), 
cone(nj), size_min(n_i), size_max(n_i), stochas- 
tic(nj), ((xJnj(nJ,nx), y jnj(nj,nx), zjnj(n j,nx), 
Z-inj(nJ,nx), flowf(n-i,nx), v_inj(n j,nx), 

alpha_inj(n j,nx), beta_inj(n j,nx), theta _inj(n j,nx), 
dtheta Jnj (n j,nx) , s wlr .angle (n _i, nx) ) , 

nx= 1 ,no_of_holes(n_i) ) , (atom Type (n J,nx) , 

breakup _type(n_i,nx), diaJrole(n j,nx), 

delp_inj(nj,nx), liq-vel(n j,nx), gas_u(n j,nx), 

gas_v(n j,nx), gas_w(n j,nx), pel-start (n J,nx), 

pcl_end(n j,nx) ) , nx=l,no_of_holes(n j)) . 

The initial droplet distribution for a given in- 
jector could be specified by making use of one of the 


In this option, we need to specify several pa- 
rameters including no_of_holes(), no_of_streams(), & 
no_of_droplet_groups(). And depending on what is 
specified for the input of the integer variable, lmdis, in 
ncc jnjector.in.l, the droplet-size distribution within 
each one of the streams is determined based on the 
following three choices: 

• In one choice, it is calculated based on a corre- 
lation typical of those widely used in describing 
the initial droplet size distribution [4]: 


dn 

n 


= 4.21 x 10 6 


_d_ 

<^32 


3.5 


r 16 - 98 ^) 0 ' 4 — 

ds2 


(67) 


where n is the total number of droplets and dn is 
the number of droplets in the size range between 
d and d + dd. This correlation also requires the 
specification of Sauter mean diameter, d.32 . Fig. 
5 shows the droplet size distribution generated 
by Eq. (67) for a case studied in [10]. The solid 
line shows the droplet number variation versus 
drop size and the dashed line shows the inte- 
grated mass variation with drop size. The drop 
size distribution within the spray is represented 
by a finite number of droplet classes as given by 
the variable, no_of_droplet_groups(). 
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Table 6. Input file: ncc_liquidjnjector.in.01. 

Input file content 

comments 

heading 

title of a description of this file 

heading 

title of property 

nolc(nJ) 

denotes the total number of initial liquid components 
in the nJ-th injector. 

heading 

title of property 

out _string 

chemical symbols of initial liquid species: e.g., C6H12. 

(ymki(n_i,n_l), 

nJ=l,nolc(ni)) 

initial mass fractions of liquid species 

heading 

title of controlling parameters 

tdrop(n_i), 
isuperhd(n J) 

tdrop(nJ) denotes the initial droplet temperature. 

If isuperhd(n J) = .true., it invokes the superheat vaporization 
model. Otherwise it makes use of a normal vaporization model. 

heading 

title of controlling parameters 

atomization(n J) , 
drop_breakup_ 
model(nJ), 
spray_table(n_i), 
steady _spray unodel 

If atomization(n_i) = .true., the initial droplet conditions are 
determined from the use of a primary atomization model. 

Otherwise they are determined from either known 
experimental conditions or a widely-used correlation. 

If drop-breakup_model(n j) = .true., it invokes the secondary 
droplet breakup option. 

If spray -table = .true., initial droplet location, velocity, size, 
and mass flow rate are input through the file - nccJiquid-table.in.l. 
Otherwise they are determined based on some form of specified 
spray correlations and configurations. 

If steady _spray -model = .true., it invokes a steady state 
spray model commonly used in many spray codes, where 
whenever the spray solver is called, after first introducing 
a new group of spray particles, it continues with the 
liquid-phase computation until after all the particles are taken 
out of the computational domain (note: It is NOT 
recommended to use this option as a steady state calculation 
could be better arrived at by making use of some features of 
the unsteady option). The solution from the unsteady option 
(steady_spray .model = .false.) is determined based on the 
values assigned to the controlling time steps - dtml, dtil, 

& dtgl, which are internal to the spray solver. 

heading 

title of controlling parameters 

no_of_holes (n_i) , 
no_of_streams(nj), & 
no _of _drople t -groups (n J ) 

no_oLholes(nJ) = The number of holes per injector (note: 
when atomizationQ = .true, or spray -table = .true., the 
no_of_holes(nj) is set equal to 1). 
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Table 6. Input file: ncc_liquidJnjector.in.01 (continued). 


no_of_streams(n J) = The number of droplet streams per hole. This 
variable is introduced mainly to distinguish the variation in the 
droplet groups based on angular orientation (note: when atomizationQ 
.true, or spray .table = .true., the no_of_streams is set equal to 1). 

no_of_droplet_groups(n j) = The number of droplet groups per 
stream. This variable is introduced mainly to distinguish the 
droplet groups based on the droplet-size variation (note: 
when spray .table = .true, or atomizationQ = .true., all different 
injected droplets are lumped together into a single variable, 
no_of_droplet_groups()) . 

heading 

title of controlling parameters 

lmdis(nJ), 
smdm(nJ), 
cone(nJ), 
size_min(n j), 
size_max(n_i), & 
stochastic(n_i) 

If lmdis(nJ) = 1, it invokes the droplet size distribution as given 
by Eq. (67); If lmdis(nJ) = 2, it assumes a uniform droplet distribution 
between the maximum and minimum limits as specified by size_max() 

& sizejninQ; & If lmdis(nJ) = 3, it invokes the droplet size 
distribution as given by Eq. (68). 

smdm(nJ) = Sauter mean diameter 

If cone(nJ) = .true., it activates a 3D solid or hollow cone spray 
configuration as shown in Fig. 2. Otherwise it activates a 2D configuration 
depending on the value chosen for the logical variable - axisymmetric. If 
it is set equal to .true., it invokes the axis-of-symmetry case as shown 
in Fig. 3 otherwise it invokes the planar case as in Fig. 4. 

size_min(n_i) & size_max(nJ) = The variables are associated with 
the lmdis(nJ) = 2 option. 

stochastic (n_i) is used in conjunction with Option (2) of Sec. 9. If set 
equal to .true., it introduces some randomness into the determination 
of initial droplet velocity distribution. 

heading 

title of controlling parameters 

((x_inj(n_i,nx), 

yinj(n_i,nx), 

z_inj(nJ,nx), 

flowf(n J,nx), 

vJnj(n_i,nx), 

alpha_inj (n_i,nx) , 

beta_inj(nJ,nx), 

(xjnj(n_i,nx),y Jnj(n_i,nx),z_inj(nJ,nx)) = spatial coordinates 
of the initial injector location. 

flowf(n_i,nx) = injector mass flow rate, (units - kgm/s for 3D & 
axis-of-symmetry & kgm/s/m for 2D planar. REMEMBER FOR AXIS-OF- 
SYMMETRY, IT IS THE TOTAL FLOW RATE OVER 360 DEGS.) 
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Table 6. Input file: nccJiquid_injector.in.01 (continued). 

theta jnj (nj,nx), 
dtheta jnj (n j,nx) , 
swlr_angle(n j,nx) ) , 
nx=l,no_oOioles(n j)) 

v_inj(nJ,nx) = droplet injection velocity, m/s. 
alpha Jnj (n_i,nx) = angle of rotation from the x-y plane, 
beta Jnj (nj,nx) = angle of rotation from the x-axis. 
theta jnj (nj,nx) = cone angle. 

dtheta jnj (n j,nx) = half-cone angle (note: Although dtheta jnj (n j,nx) 
= theta jnj (n j,nx)/2 for SOLID CONE SPRAY, it is invoked by 
setting dtheta jnj (nj,nx) either equal to 0 or theta jnj (nj,nx)/2) . Refer 
to Figs. 2-4, for a better understanding of the angular representation. 

swlr_angle(n j,nx) = swirl angle. 

heading 

title of controlling parameters 

(atom_type(n j,nx), 
breakup _type(n_i,nx), 
dia holefn i,nx), 
delp jnj (n_i,nx) , 
liq_vel(n j,nx), 
gas_u(n j,nx), 
gas_v(nJ,nx), 
gas_w(n j,nx), 
pcLstart(nJ,nx), 
pcl_end(nJ,nx)), 
nx= 1 , no _of Jroles (n j ) ) 

atom jype(n j,nx) = 0 => no primary atomization specified, 

= 1 => blob jet primary atomization model, 

= 2 => sheet breakup primary atomization model, 

= 3 => air blast primary atomization model, 

= 4 => BLS primary atomization model. 

breakup jype(n j,nx) = 0 => no secondary droplet breakup model specified, 
= 1 => Rayleigh-Taylor secondary droplet breakup model, 

= 2 => TAB secondary droplet breakup model, 

= 3 => ETAB secondary droplet breakup model. 

diaJrole(n j,nx) represents the injector orifice diameter (m). 

delp jnj (nj,nx) represents the pressure drop across the injector (Pa) 

(note: needed for pressure swirl atomization only). 

liq_vel(n j,nx) represents the liquid velocity from annulus (m/s) 

(note: needed for air-blast atomization only). 

gas_u(n j,nx), gas_v(n j,nx), gas_w(n j,nx): gas velocity components 
in annulus adjacent to liquid film (m/s) (note: needed for air-blast 
atomization only). 

pcl_start(n j,nx): starting parcel number for the annular air-blast injector. 
pcl_end(n j,nx)): ending parcel number for the annular air-blast injector. 
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Fig 2a. Geometrical details of fuel injection for 
a 3D solid or hollow cone spray. 



§ = azimutal angle, deg 
number of rays = 8 ((|)=45) 
For a 3D cone, 
no_o f_st reams ( ) should 
defined as a multiple of 8 
number along each ray = 
no_of_st reams 0/8. 


Fig 2b. Initial spray particle orientation in a 
circular cross-section. 
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9 = cone angle, deg 
d0= half-cone angle 



Fig 3a. Geometrical details of fuel injection for 
an axis-of-symmetry case. 


• 00 

* 

• 01 


For axis of symmetry: 
no_of_streams ( ) is distributed 
between 01 to 00 . 


; C 


Fig 3b. Initial spray particle concentration in an 
axis-of-symmetry case. 
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Table 7. Input file: nccj3pray_table.hr. 01. 

Input file content 

comments 

nos(n_i) 

denotes the total number of droplet groups in the n_i-th 
injector. 

(ni,xxJnj,yyJnj, 

ni = number of the droplet group. Its value ranges between 

zz _inj ,uu Jnj , vv_inj , 
ww _inj ,r_inj ,fld_d) 

1 to nos(nJ). 

,ni=l,nos(nJ)) 

(xx Jnj,yy_inj,zz_inj) = spatial coordinates. 
(uu_inj,vvJnj,ww_inj) = velocity components. 
r_inj = droplet size in radius. 

flcLd = mass flow rate of the ni-th droplet group (note: 
SUMMATION OF fld_d OVER ALL THE nos(n_i) 
DROPLET GROUPS IS EQUAL TO THE TOTAL MASS 
FLOW RATE OF THE INJECTOR). 



Figure 5 Droplet-size distribution. 


• In the second choice, the initial droplet sizes are 
distributed evenly between two specified maxi- 
mum and minimum drop sizes - size_max() & 
size rriiiiQ as specified in Table 6. And the size 
interval depends on the value specified for the 
variable, no_of_droplet_groups(). However, this 
option is applicable only in certain special cases. 

• In the third choice, the initial droplet size dis- 
tribution is calculated by making use of a cliped 
probability distribution function as given by 

P(X) = [ f( x )dx/ [ f(x)dx (68) 

Jo Jo 

where f(x) = (1.0 — exp(— x){l+x+\x 2 + \x ? '). 
After a random sampling of 0 < % < 12.0, the 
initial droplet diameter is determined as follows: 
(0.3 < 4 P(x) < 1.5)d 3 2. If the value for 4P(y) 
goes beyond the limits of 0.3 and 1.5, the ran- 
dom sampling is reiterated until that number 
falls within the range of the lower and upper 
bounds. 

Depending on what is specified for the logical 
variable, cone, of Table 6, the droplet velocity dis- 
tribution amongst various streams of a given hole is 
calculated by assuming the spray to be either a solid 
or hollow cone spray. A graphical illustration of three 
different cone configurations are shown in Figs. 2 to 
4. Figs. 2a & 2b refer to a 3D case, Figs. 3a & 3b 


NASA/CR— 2008-215290 


21 





to an axisymmetric case, and Figs. 4a & 4b to a 2D 
planar case. It is noteworthy that in an axisymmet- 
ric case, the x-axis is assumed to be aligned with the 
axis-of-symmetry. 

Fig. 2a shows the geometric details of a hollow 
cone spray in 3D where 0 is known as the cone angle, 
dO is the half-cone angle, and for a solid cone spray 
dO = d/2. And a represents the angular rotation 
from the x-y plane and (3 is the angular rotation from 
the x-axis. Fig. 2b shows initial spray stream orien- 
tation in a circular cross section. We try to simulate 
the spray by a finite number of streams as given by 
the variable, no_of_streams(). Each one of the dark 
circles in Fig. 2b represents a different stream. The 
streams are distributed evenly along each one of the 
different rays which are separated from each other 
with an angle of separation as given by <f> . Currently, 
we have hard-coded the angle of separation to be 45° 
which means we have restricted the number of rays 
to be eight in 3D. Therefore, when specifying a num- 
ber for no_of_streams(), it should be borne in mind 
that it should be a multiple of eight. And the num- 
ber of streams along each one of the rays, therefore, 
becomes equal to no_of_streams()/8. In Fig. 2b, the 
no_of_streams() has a value of 32 and, therefore, the 
number of streams along each one of the rays for this 
case is 4. 

Fig. 3a. shows the geometric details for an 
axisymmetric case. Since the computations are per- 
formed only in the first quadrant of the x-y plane, all 
of the specified number of streams, no_of_streams(), 
are distributed evenly between the lower and upper 
halfs as shown in Fig. 3a and 3b. Here, the upper half 
refers to the region between 01 to 00 and the lower 
half refers to LO to LI. It is noteworthy that each 
droplet group in a given stream represents a circular 
ring of liquid for an axisymmetric case. 

The geometric details for a 2D case are shown 
in Figs. 4a & 4b. Here, it is noteworthy that 
each droplet group in a given stream represents a 
planar sheet of liquid. All the specified number of 
streams, no_of_streams(), are distributed evenly over 
both sides of the cone center as shown in Fig. 4b. 

For each one of the different injector holes, 
no_of_holes(), it requires the specification of the fol- 
lowing parameters as described in nccJnjector.in.l 
of Table 6 (However, because of the geometric dif- 
ferences that exist between 3D, 2D, and axisymmet- 
ric configurations, some of the input parameters may 
have different units as noted below): 


• The initial (x, y, z) coordinates of the hole loca- 
tion. 

• The mass flow rate per hole - however, the defi- 
nition of the units for the injector mass flow rate 
per hole differs: it is kgm/s for 3D & the axis- 
of-symmetry and kgm/s/m for 2D planar (note: 
In an axis-of-symmetry, the specified mass flow 
rate per hole refers to the entire mass flow rate 
over 360 degrees). 

• The following variables define the angular ori- 
entation: cti n j =angle of rotation from the x-y 
plane, /3i n j = angle of rotation from the x-axis, 
0 in j = cone angle, & 00j nj - = half-cone angle 
(note: Although dd in j = 0 inj -/ 2 for a solid cone 
spray, a specified value of zero for d9i n j also in- 
vokes a solid-cone spray configuration). 

• The variable, swlr_angle(), allows a means to 
specify the tangential component in the case of 
both 3D and axisymmetric sprays. 

Stochastic Option 

In the above description of Option (2), the ini- 
tial droplet conditions are determined based on a de- 
terministic approach. However, such an approach 
when employed in 3D may lead to time consuming 
calculations requiring the use of many droplet groups. 
With this in mind, we wanted to introduce some ran- 
domness into the calculation of the initial droplet ve- 
locity distribution as an option. This option can be 
invoked by setting the logical variable, stochastic (), 
to be .true. . In a stochastic approach, the initial an- 
gular orientation of a particle could be randomized in 
a number of different ways. The obvious way would 
be to determine the initial droplet direction based on 
a complete randomization in a field as defined by the 
geometric parameters: solid cone angle, 0, half cone 
angle, dO and/or azimuthal angle, <j>. Another way 
is try to retain the basic features of the determinis- 
tic approach where the initial droplet directions are 
taken as the means in a random field with the ran- 
domness localized to within the surrounding angular 
sub-segment of the respective droplet as represented 
by 66, 6dd, SdScj), or 5d95(j>. We have chosen to go 
with the later option as it preserves the same basic 
structure but provides some benefits of the stochastic 
approach. 
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Option (3): 


Here many aspects of fuel injection remain the 
same as in Option (2) but the initial droplet size 
and velocity distribution may differ depending upon 
the primary atomization model selected. More de- 
tails on the available primary atomization models 
can be found in the appendix. For each one of the 
different injector holes, no_of_holes() it requires the 
specification of all the parameters as in Option (2). 
However, it doesn’t make use of the values set for 
size_max() & size_min(). Also, it requires the spec- 
ification of some additional parameters as described 
in ncc_injector.in.l of Table 6. However, there exist 
some differences with regard to their usage depending 
upon the primary atomization model. Some major 
differences are highlighted below: 

• The integer variable, atom_type() of Table 6, is 
used to select the primary atomization model of 
one’s choice. 

• dia_hole() needs to be specified for all primary 
atomization models. 

• delp_inj() needs to be specified only with the 
sheet breakup primary atomization model and 
it is used in Eq. (9) of the appendix. 

• liq_vel(nJ,nx), gas_u(nJ,nx), 

gas_v(nJ,nx), gas_w(n _i,nx), pcl_start(n J,nx), & 
pcl_end(n j,nx)) are required only for the air- 
blast primary atomization model. 

• In all of the primary atomization models except 
for the air-blast, the total number of droplet 
groups (particles) to be injected at the time 
of every injection period is specified though 
the use of the variable, no_of_clroplet_groups(). 
For the air-blast, it is given by pcl_end(n_i,nx)- 
pcl_start(n j,nx)+l. 

• In this option, all of the injected droplet groups 
as specified by either no _of -droplet -groups () 
or pcLend(n-i,nx)-pcl_start(n_i,nx)-|-l are dis- 
tributed randomly across the entire cross-section 
of any selected 2D, axisymmetric, or 3D spray 
configuration. So it is recommended that 
no_of_streams() be set equal to one as no dis- 
tinction is made between streams. 

The secondary droplet breakup option can be 
invoked regardless of how the initial droplet condi- 
tions are generated using Options (1) to (3). The sec- 
ondary droplet breakup option can be invoked by set- 
ting the logical variable, drop -breakup _model(), to be 



Fig. 6 A vector illustration used in the particle search 
analysis. 

.true.. And then the integer variable, breakup_type() 
of Table 6, is used to select the secondary droplet 
breakup model of one’s choice. More details on the 
available secondary droplet breakup models can be 
found in the appendix. 

10 SPRAY SOLUTION ALGORITHM 

In order to evaluate the initial conditions that 
are needed in the integration of the liquid-phase equa- 
tions, we first need to know the surrounding gas- 
phase properties at each particle location. But in or- 
der to evaluate the gas-phase properties, it is first nec- 
essary to identify the computational cell in which a 
given particle is located. It is a trivial task to track a 
particle in the regular rectangular coordinates. How- 
ever, the particle tracking becomes complicated when 
the computational cell is no longer rectangular in the 
physical domain and it becomes even more compli- 
cated when the particle search is performed within 
the context of parallel computing. 

We have developed and implemented an effi- 
cient particle-tracking algorithm for use with parallel 
computing in an unstructured grid. The search is ini- 
tiated in the form of a local search originating from 
the computational cell in which the same particle was 
found to be located in the previous time step. The 
location of the computational cell is then determined 
by first evaluating the dot product of x pc . a n = x pc 
\a n \ cos (</>), where x pc is the vector defined by the 
distance between the particle location and the center 
of the n-face of the computational cell, a n is the out- 
ward area normal of the n-face as shown in Fig. 6, 
and <f> is the angle between the two vectors. 

A simple test for the particle location requires 
that the dot product be negative over each and ev- 
ery one of the n-faces of the computational cell. If 
the test fails, the particle search is carried over to the 
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adjacent cells of those faces for which the dot prod- 
uct turns out to be positive. Some of those n-faces 
might represent the boundaries of the computational 
domain while the others represent the interfaces be- 
tween two adjoining interior cells. The search is first 
carried over to the adjacent interior cells in the di- 
rection pointed out by the positive sign of the dot 
products. The boundary conditions are only imple- 
mented after making sure that all other remaining 
possibilities point towards a search exterior of the 
computational domain. This implementation ensures 
against any inadvertent application of the boundary 
conditions before the correct interior cell could be 
identified. 

After the gas-phase properties at the particle 
location are known, the solution for the ordinary dif- 
ferential equations of particle position, size, and ve- 
locity are advanced by making use of a second-order 
accurate Runge-Kutta method. The partial differen- 
tial equations governing the droplet internal thermal 
and mass transport are integrated by making use of 
a fully implicit Newton-Raphson iteration method. 

Finally, the liquid-phase source terms of the gas- 
phase conservation equations (1-4) are evaluated by 
making use of a time-averaging method. 

11 DESCRIPTION OF COMPUTER 
FLOWCHART, & THE TIME- AVERAGING 
SCHEME USED IN THE GAS-PHASE 

SOURCE-TERM CALCULATIONS 

• In order to know more about the time-averaging 
method, we need to know first about the three 
different time steps that are internal to the spray 
code. Af m /, Atui and At g i. 

A t m i - the actual time step used in integrating 
the liquid-phase equations which is determined 
based on the smallest of the different time scales 
associated with various rate-controlling phenom- 
ena of a rapidly vaporizing droplet. Some of 
the limiting scales are the average droplet life- 
time, the time it takes for the droplet to traverse 
the local grid spacing, the time it takes for the 
droplet internal temperature to reach the liquid 
boiling temperature, & a relaxation time scale 
associated with droplet drag. The time-scale re- 
striction based on these criteria can become quite 
severe for smaller drops - for drops of sizes less 
than a micron. 

Atu - the injection time step. It is the time 
step at which a new discretized parcel of different 


droplet groups are introduced into the computa- 
tion. 

A t g i - the global time step. Its introduction 
seems to provide better convergence in both un- 
steady and steady-state computations. 

• When the spray solver is called it advances the 
liquid phase equations over a number of itera- 
tions as determined by the ratio of At g i/ At m i- 

• It then evaluates the time-averaged contribution 
of the liquid-phase source terms, S g i , of the gas- 
phase governing equations (1-4) as follows: 

- E ^ *- <«> 

m = 1 y 

where 


M 

= A t gl (70) 

771— 1 

• The values for A t m i, A t g i, & Atu are specified 
in the input file, nccJiquid_solver.in, of Table 8. 

• In steady-state computations, it is recommended 
to use for both A t g i and Atu a value of about 
1 ms which is roughly equivalent to the aver- 
age lifetime of the droplets for a typical reacting 
spray encountered in conventional low-pressure 
gas-turbine combustors. 

The averaging scheme could be explained better 
through the use of a flow chart shown in Fig. 7. 
The main spray solver is invoked with a controlling 
routine, DCLR, which, then, executes the following 
steps: 

1. It first initializes the source terms to zero. 

2. Updates the global time, t g i, based on A t g i. 

3. Checks to see if t m i < t g i < t m i + A t m i. If it 
is, it returns control over to the calling routine 
and supplies the other flow solvers, e.g., flow or 
EUPDF, with the source terms, S g i, of Eqs. (1)- 
(4). If not, it proceeds with the next step. 

4. Checks to see if it it is time to introduce a new 
group of particles. 
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Fig. 7 The flowchart of LSPRAY-lll. 
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Table 8. Input file: nccJiquid_solver.in. 

Input file content 

comments 

heading 

title of controlling parameters 

ldread,ispray_mod, 

( when_st art _spr ay (n j ) , 
n J=l,no_of_injectors) 

If ldread = .true., restarts the calculation from the data stored from 
the previous computation. Otherwise initiates a new spray computation. 

ispray_mod controls the calls to the spray solver. The spray solver is called 
once at every other CFD iteration as given by the number, ispraycmod. 

when_start_spray() represents the CFD iteration where the computations for 
the n_ith injector are initiated. 

heading 

title of controlling parameters 

dtml, dtgl, & 
dtil 

dtml = time step for advancing the liquid phase equations. 

dtgl = global time step. Whenever the spray solver is called, it advances 
the spray computations over a period of dtgl before returning control 
over to the calling routine. To be more precise, it advances the 
liquid phase equations over a number of time steps as determined by 
(dtgl/dtml). 

dtil = injection time step. It is where a new group of droplets are 
introduced into the computation. 


5. Proceeds with solving the liquid-phase equations 
with calls to the following routines: 

• Calls the particle tracking routine and as- 
signs particles based on the parallel strategy 
implemented. 

• Interpolates gas-phase properties at the 
particle location. 

• Advances liquid-phase equations and, then, 
deletes any particles that are no longer 
needed in the computations. 

6. Evaluates the liquid-phase source-term contribu- 
tions, S m i, of Eq. 69. 

7. Updates the time, t rn [, based on A t m i. 

8. It then goes back to step (3) and repeats 
the whole process again until the computa- 
tions are completed over a global time step: 

Ern=l = Atgl. 

12 IMPLEMENTATION OF BOUNDARY 
CONDITIONS 


• In case of the droplet impingement with the com- 
bustor walls, the droplets may evaporate, move 
along the wall surfaces, and/or reflect with re- 
duced momentum. The physics of droplet im- 
pingement with the walls is not completely un- 
derstood. In our present calculations, it is as- 
sumed that the droplets, after having lost most 
of their momentum upon impingement with the 

* walls, move along the wall surfaces with a veloc- 
ity equal to that of the surrounding gas. 

The implementation of the periodic boundary 
conditions becomes rather complicated as it re- 
quires taking into consideration several aspects 
arising from a particle leaving the computational 
domain from one periodic boundary has to reen- 
ter the domain back from a corresponding sec- 
ond periodic boundary with appropriate condi- 
tions. The periodic boundary conditions are im- 
plemented with the help of some appropriately 
defined transformation matrices. 


The spray code supports most of the boundary 
conditions that are in use in the current version of 
the NCC CFD module but not all. Here, we would 
like to highlight some of the boundary conditions: 
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• The symmetric boundary condition is imple- 
mented in such a way to satisfy the criterion that 
for every particle crosssing the symmetry line, a 
similar one re-enters the domain in a direction 
given by the reflection off of the symmetry line. 

• When the particles try to move out of the exit 
boundary, they are removed out of the compu- 
tation. 


13 DETAILS OF COUPLING BETWEEN 
LSPRAY-III AND OTHER SOLVERS 

• The spray code is designed to be a stand-alone 
module such that it could easily be coupled with 
any other unstructured-grid CFD solver and the 
same holds true for EUPDF. (However, some 
grid-related parameters such as area vectors, grid 
connectivity parameters, etc. need to be sup- 
plied separately). 

• The spray solver needs information on the gas 
velocity and scalar fields from the other solvers 
and, then, it in turn supplies the liquid-phase 
source terms. 

• The PDF solver needs information on the mean 
gas velocity, turbulent diffusivity and frequency 
from the CFD solver and the liquid-phase source 
terms from the spray solver, and then it in turn 
provides the solution for the scalar (species and 
energy) fields to the flow and spray solvers. 

• It should also be noted that both the PDF and 
spray solvers are called once at every other spec- 
ified number of CFD iterations. 

• All of the three solvers (LSPRAY-III, EUPDF, 
and CFD) are advanced sequentially in an iter- 
ative manner until a converged solution is ob- 
tained. 

• All three codes (EUPDF, CFD, and LSPRAY- 
III) were coupled and parallelized in such a way 
to achieve maximum efficiency. 

The coupling issues could be better understood 
through the use of a flow chart shown in Fig. 8. 
It shows the overall flow structure of the combined 
CFD, LSPRAY-III, and EUPDF modules. Both the 
PDF and spray codes are loosely coupled with the 


CFD code. The spray code is designed in such a 
way that only a minimal amount of effort is needed 
for its coupling with the flow and PDF solvers. The 
present version of the spray module relies entirely on 
the use of Fortran common blocks for its informa- 
tion exchange with other modules. Even this reliance 
should entail only few changes to be made within the 
spray code for its linkage with different solvers. The 
PDF code is also structured along similar coupling 
principles. 

The flow chart of Fig. 8 contains several blocks 
- some shown in black and/or solid lines and the oth- 
ers in color and/or dashed lines. The ones in solid 
blocks represent the flow chart that is typical of a 
CFD solver. The ones in dashed blocks represent the 
additions arising from the coupling of the spray and 
PDF solvers. 

The coupling starts with the calling of the 
subroutine - spray _pdf_read_parameters, which 
then reads the spray control parameters from the 
input file, nccdiquidsolver.in of Table 8. This 
table provides a detailed description of the fol- 
lowing input file variables - ldread, ispray_mod, 
(when_start_spray(n_i), nJ=l,no_of_injectors), dtml, 
dtgl, & dtil. The coupling is then followed by the 
calling of the pdf int rerun subroutine. It initial- 
izes PDF computations and, also, it may restart 
the PDF computations if needed from the data 
stored from a previous iteration. Similarly, we 
call spray _int_rerun for the spray computations. 
It is noteworthy that the spray computations can 
be restarted by reading the data from the restart 
files: nccJiquid-params.out , ncc_liquid_results.db , 
& nccJtiquid_results.mc.db. The restart capability 
is invoked by setting the logical variable, ldread, of 
the input file, nccJliquid.solver.in of Table 8, to be 
true. Otherwise, the spray computations are initial- 
ized to start from the beginning. Then, the coupling 
proceeds with the calling of the following subroutines: 
dclr for integrating the spray calculations and eupdf 
for the Monte Carlo PDF. The input variable, is- 
prayunod of Table 8, controls the calls to the spray 
integrating routine. The spray solver is called once at 
every other number of CFD iterations as specified by 
ispray_mod. And the first call to the spray solver is 
controlled by the input variable, when_start_spray() 
of Table 8, which represents the starting CFD itera- 
tion number from which the spray computations are 
initiated. Finally, the coupling ends with the call- 
ing of a subroutine, spray _pdf_output, which will 
create a set of new restart files. 
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START 


SPAWN MULTIPLE PEs FOR PARALLEL COMPUTING 

i 

READ PARAMETERS & GEOMETRIC DATA 

, -L- 

I CALL SPRAY PDF READ PARAMETERS I 


INITIALIZE 

InitialIze^pdf VARIABLES | 


c CALL PDF INT RERUN ^ 




^[rerun 


READ PDF RESTART FILES 


zJ- 



Fig. 8 The overall flowchart of the combined CFD, LSPRAY-III, and EUPDF solvers. 
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14 DETAILS OF FORTRAN 
SUBROUTINES & FUNCTIONS 

Table 9 provides a list of all the Fortran sub- 
routines developed as a part of the spray module. 
This table also provides information on all the For- 
tran functions. It also describes the purpose of all 
the individual subroutines and functions. Similary, 
Table 10 provides a list of all the Fortran functions 
and subroutines associated with the implementation 
of the high pressure equation of state as described in 
Sec. 3, and the high pressure behavior on the gas- 
phase transport properties as in Sec. 4. 

15 DETAILS OF PARALLELIZATION 

There are several issues associated with the par- 
allelization of both the spray & PDF computations. 
The goal of the parallel implementation is to extract 
maximum parallelism so as to minimize the execu- 
tion time for a given application on a specified num- 
ber of processors [37]. Several types of overhead costs 
are associated with parallel implementation which in- 
clude data dependency, communication, load imbal- 
ance, arithmetic, and memory overheads. The term 
arithmetic overhead is the extra arithmetic opera- 
tions required by the parallel implementation. Mem- 
ory overhead refers to the extra memory needed. Ex- 
cessive memory overhead reduces the size of a prob- 
lem that can be run on a given system and the 
other overheads result in performance degradation 
[37]. Any given application usually consists of sev- 
eral different phases that must be performed in cer- 
tain sequential order. The degree of parallelism and 
data dependencies associated with each of the sub- 
tasks can vary widely [37]. The goal is to achieve 
maximum efficiency with a reasonable programming 
effort [37]. 

In our earlier work, we discussed the parallel 
implementation of a spray algorithm developed for 
the structured grid calculations on a Cray T3D [10]. 
These computations were performed in conjunction 
the Monte Carlo PDF method. The parallel algo- 
rithm made use of the shared memory constructs ex- 
clusive to Cray MPP (Massively Parallel Processing) 
Fortran and the computations showed a reasonable 
degree of parallel performance when they were per- 
formed on a NASA LeRC Cray T3D with the number 
of processors ranging between 8 to 32 [10]. Later on, 
the extension of this method to unstructured grids 
and parallel computing written in Fortran 77 with 


PVM or MPI calls was reported in [11-15]. The lat- 
est version in Fortran 77/90 offers greater computer 
platform independence. In this section, we only high- 
light some important aspects of parallelization from 
Refs. [10-15], 

• Both the EUPDF and CFD modules are well 
suited for parallel implementation. For the gas- 
phase computations, the domain of computation 
is simply divided into n-Parts of nearly equal 
size and each part is solved by a different pro- 
cessor. Fig. 9a illustrates a simple example of 
the domain decomposition strategy adopted for 
the gas-phase computations where the total do- 
main is simply divided equally amongst the avail- 
able computer processing elements (PEs). In 
this case, we assumed the number of available 
PEs to be equal to four. 

• But the spray computations are more difficult to 
parallelize for the reasons summarized below: 

(1) Non-uniform nature of spray distribution: 
Most of the particles are usually confined to a 
small region where the atomizer is located. 

(2) Dynamic nature of Lagrangian particles: 
Particles keep moving between different sub- 
divided domains of an Eulerian grid (grid used 
in the CFD computation) which are assigned to 
different processors. While some new particles 
are introduced at the time of fuel injection, some 
others are taken out of computation. 

Conceptually, there are several ways to paral- 
lelize the spray computations, we, however, devel- 
oped and tested two different domain decomposition 
strategies [10-14]. 

• Strategy I: 

The Lagrangian particles were assigned fairly 
uniform amongst the available processors but the 
calculations associated with the particle track- 
ing, the interpolation of the gas-phase proper- 
ties, and the source-term evaluation were com- 
puted on the processor of the computational grid 
in which a particle is located. 

This strategy leads to an uniform loading during 
integration but leads to excessive message pass- 
ing. 

• Strategy II: 
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Table 9. Description of LSPRAY-III Fortran subroutines & functions. 

Function 

Purpose of the Function 

blasiu(x) 

This function returns a solution for the function, f(Bk), 
of Eq. (21) for use in computing the droplet regression rate. 

Subroutine 

Purpose of the Subroutines 

axisym_part _adj ust 

This routine was written by Jeff Moder and essentially 
performs the same functions as the intpla routine but 
modified for application with a 2d-axisymmetric case. 

calc_F_suvw() 

It evaluates the right-hand sides of ODEs for both 
droplet size and velocities, and then provides the 
solution to the calling routine, chaslv. 

chaslv 

This routine has two main functions: 

(1) It integrates the liquid phase equations. 

(2) It removes the particles that are no longer needed 
in the computation. 

dclr 

This routine is called once at every other CFD iteration as 
specified by ispray_mod. It is primarily a controlling 
routine for spray computations. 

This routine has the following functions: 

(1) It initializes the source terms to zero. 

(2) Checks to see if new particles need to be introduced. 

(3) Advances liquid phase equations over an allowable or 
pre-specified time step, dtml, with calls to the following 
routines: 

intplal - Interpolates the gas phase properties at the particle 
location. 

chaslv - Advances liquid phase equations. 

intpla - Identifies computational cells and PEs associated 
with spray particles. 

sprips - Evaluates the liquid phase source term contributions 
of the CFD and PDF equations. 

(4) (a) For unsteady spray computations (steady_spray_model=F), 
this routine is called once at the beginning of every 

global time-step, dtgl. This is accomplished by continuing 
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Table 9. Description of LSPRAY-III Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 


with steps (2) and (3) until the computations are completed 
over one single global time step, cltgl. (b) For steady spray 
spray computations (steady .spray _model=T), this routine is 
only called once to solve entire lifetime of all particles 
injected. It repeats step (3) until there are no more spray 
particles left to be integrated (note: The particles are taken 
of computation when they reach either a certain size of 
negligible proportion or exit out of the computational domain). 

(5) Returns control over to other solvers, e.g. CFD and EUPDF, 
and supply them with the source terms averaged over dtgl. 

dropdis(rhol, 

flowdum,sr, 

fld,smd,nofg) 

This routine provides initial droplet distribution from 
the following correlation (also, appears as Eq. 67): 

dn/n = a((D/D 32 ) al P)exp(-b((D/D 32 ) bet ))dD/D 32 
where a, b, alp, and bet are constants. 

drop jc(n_i, 
nmih,nmis, 
nmip) 

It is called by the spray jnainJnjection routine in conjunction 
with spray _table(n_i))=F at the beginning of every injection 
time step, dtil, for introducing a new group of spray particles. 

It applies for the introduction of new particles in conjunction with 
Options (2) or (3) of Sec. 9. In Opt. (2) the initial conditions 
are specified from the known correlations and in Opt. (3) from 
the available primary atomization models. 

find_transport _ds ( 
ijle,chem_tnodel, 
number _of_species,hp_eos, 
temp gas, pres gas.y gas, 
rhogm,visgm,congm,difgm) 

It computes the following properties of a gas mixture 
at the droplet interface by making use of the one-third 
rule of Eq. (46): molecular viscosity, gas density, 
thermal conductivity, and diffusion coefficient. 

find_xyzface 

(i,xcfac,ycfac, 

zcfac) 

This routine computes x, y, and z locations of all the face 
centers of an element, i. This information is used in the 
particle search algorithm. 

get _liq_t at .properties 
(tmp,pmp,j) 

It computes the following variable properties of a liquid 
mixture: density, specific heat, molecular viscosity, 
gas density, thermal conductivity, and diffusion coefficient. 
This routine also computes the surface tension of a 
multi-component liquid fuel. 

get _liq_t at .properties 1 

(tmp,pmp,boil_tem, 

cp_sph,hvap_sph,icuin) 

This subroutine is used in conjunction with the superheat 
vaporization model and computes the following variable 
properties of the liquid fuel: boiling temperature, specific 
heat at constant pressure, surface tension, heat of vaporiza- 
tion, molecular viscosity, & liquid density and volume. 

get_liq_tat_properties_sp 

(tmpl,tmp2,cpltl, 

cplt2,j) 

This subroutine is used in conjunction with the superheat 
vaporization model and computes the following variable 
properties of the liquid fuel: specific heat at droplet 
temperature, & specific heat at boiling temperature. 
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Table 9. Description of LSPRAY-III Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

init_psi2_dsd 

This routine provides initial droplet distribution from 
a clipped probability distribution function as given by 
Eq. 68. 

intpla 

This routine performs three main functions: 

(1) Particle Tracking - It identifies the computational cell 
in which a particle is located. In parallel computing, 

it also means identifying the corresponding processor of the 
computational domain in which a particle is located. 

(2) It implements appropriate boundary conditions. 

(3) It reassigns the particles between different PEs based on 
the parallel strategy employed. 

intplal 

This routine interpolates the gas-phase properties at the 
particle location. In the present case, a simple first-order 
interpolation is employed. 

mimd_spray 

For computational elements where the neighboring cells are 
partitioned between different processors, this routine initializes 
the arrays, ipr_fr_id() and ile_fr_id(), by storing the information 
on the corresponding processor and element ID numbers of 
all the neighboring cells. This information is needed in order to 
track the particle movement between neighboring cells. 

mimd_spray_recv 

(i_recvfrom) 

This subroutine is called by mimd_spray in order to gather 
some relevant information from the neighboring processors. 

mimd_spray_sencl 

(i_sendto) 

This subroutine is called by mimd spray in order to send 
some relevant information to the neighboring processors. 

par_loc(xparz, 

yparz,zparz, 

ipare,iparp) 

Given the x,y,z coordinates of a particle location, the 
algorithm identifies the corresponding computational cell 
in which a particle is located. It also identifies the 
corresponding processor of the computational grid in which 
a particle is located. This information is used to locate 
newly injected droplet groups. 

prnspr 

It provides written output of the spray data in a standard format. 

spray _int_rerun 

This routine has two functions: (1) It initializes the spray 
solver based on the data read from various input files. 

(2) If it is a rerun, it restarts the computations from the 
stored data of a previous computation. 
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Table 9. Description of LSPRAY-III Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

spray .main JnjectionQ 

This routine is called to inject a new group of spray 
particles for each injector at the beginning of every 
injection time step, dtil. 

spray .output 

This routine is called to create appropriate restart files 
for the spray computations and, also, to provide for debugging 
purposes written output of some useful gas and spray data in 
a standard format. 

spray _plot_output 

It provides written output of some useful spray data in 
a standard format. It can be used in conjunction with 
steady .spray .model = F. 

spray _read_parameters 
(ipid,ncells)) 

This routine is called once in the beginning to read and 
initialize some of controlling parameters associated with 
the spray computations. 

spray .read-transport 
( (ipid, fid, filename) 

This routine is called once to read liquid thermo and 
tranport data: various constants used in the evaluation 
the polynomials in temperature for constant pressure 
specific heat, thermal conductivity, and viscosity. 

sprips 

This routine is called to provide the liquid-phase source terms of 
Eqs. (l)-(4) for use in both CFD and Monte Carlo PDF solvers: 

smlc(i) = liquid-phase contribution of Eq. (1) of Section 
2. 

smlmx(i), smlmy(i), smlmz(i) = liquid-phase contribution 
of Eq. (3) of Section 2. 

smle(i) = liquid-phase contribution of Eq. (4) of Section 2. 

sy(il,iu,bb, 

clcl,aa,cc) 

It is a tri-diagonal matrix solver. It is used in the solution 
of both Eqs. (27) & (37). 

uvw_par(swlr_angle(nx) , 
angle _work, nx, ny, nn , 
v Jnj ,nmip,t_rotation, 
cone, n_cone_r ay s , 
cone .rotation, uloc, 
vloc,wloc,n_L) 

It computes the initial particle injection velocity for different 
geometric configurations of spray as shown in Figs. 2 to 4. 
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Table 10. Description of PENG_ROB Fortran subroutines & functions. 

Function 

Purpose of the Function 

p eng _rob .density 
(chem_model,ycomp, 
tij, pressure, 
number _of_species) 

It provides a solution for the density of a multi-component 
mixture based on the Peng-Robinson equation of state. 

peng_rob .pressure 
(chem_model,ycomp, 
tij , density, 
number _of_species) 

It provides a solution for the pressure of a multi-component 
mixture based on the Peng-Robinson equation of state. 

Subroutine 

Purpose of the Subroutines 

hp.diff.table 

(ct,pt,dabcvi) 

This subroutine provides a value for {Dab P)/(D, 4 b P) + 
based on the Takahashi high-pressure correlation for the 
calculation of the diffusion coefficients in a gas mixture. 

It is based on the interpolation from the data derived from 
a graphical illustration of Reid et al [16], “The Properties 
Gases and Liquids.” The retrieved data is stored in 
hp_diff_table.dat as a table of f(tpr_diff,ppr_diff) values. 

hp .diffusion 

(tij,pij,its, 

nwarn,diin_hp) 

It is called to evaluate the gas-phase mixture-averaged 
diffusion coefficients valid at high pressures. 

hp.transport 
(calc_name,its, 
tij,pij,nwarn, 
avisi Jrp , condi _hp ) 

It is called to evaluate the gas-phase transport properties 
valid at high pressures: viscosity and conductivity. 

p eng .rob .all _rhos ( 
chem_model,ycomp, 
tij, pressure, 
number _of .species , 
rhol, rho2, rho3, 
vvl, vv2, vv3, 
zl, z2, z3, 
aprgm,bprgm) 

It provides the roots for the cubic equation, Eq. (6), involving 
the compressibility factor, Z. One of the three roots of Eq. (6) 
yields Z v and the other Z ;. Jeff Moder made significant 
improvements to this routine. 

peng_rob_gen 

(xsp,ppr,tpr, 

its,rho0) 

It takes the following as the input variables: xsp() - the mass 
fraction of the species, ppr - pressure, tpr - temperature, and 
its - the number of species. And it returns as the output, rhoO - 
density based on the solution of the Peng-Robinson EOS. 

peng.rob_zcrit.rhor 

(ppr,tpr,its) 

It provides a solution for Zj C and p lr for ith species based 
on the solution of the Peng-Robinson EOS for a given value 
of P.; r and Tj r . It also requires the specification of the 
total number of species, its. 
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Fig. 9a An illustration of the parallelization 
strategy employed in the gas flow comp- 
utations . 
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Fig. 9b An illustration of the parallelization 
strategy employed in the spray computa- 
tions . 
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Table 11. CPU time (sec) per cycle 

versus number of PEs. 



Number of processors 

Solver 

Characteristic 

2 

5 

10 

CFD 

5 steps/cycle 

2.50 

1.25 

0.75 

EUPDF 

1 step/cycle 

6.5 

2.9 

1.9 

LSPRAY 

100 steps/cycle 

1.70 

0.64 

0.53 


Table 12. CPU Time (Sec) Per Cycle Versus Number of PEs. 


Number of Processors 

Solver 

Characteristic 

1 

2 

4 

8 

16 

Spray 

100 steps/cycle 

6.83 

5.29 

2.94 

1.64 

0.87 

Max. Spray Particles in a PE 

2695 

2097 

1165 

623 

312 

Min. Spray Particles in a PE 

2695 

598 

118 

14 

0 


The Lagrangian particles were assigned to the 
processor of the computational grid where the 
particle is located. Fig. 9b illustrates a sim- 
ple example of the domain decomposition strat- 
egy adopted for the liquid-phase computations 
where the corresponding gas-flow computational 
domain is divided into equal parts between the 
four available PEs. In this strategy, the La- 
grangian particles are assigned to the processor 
of the computational grid where a particle is lo- 
cated. 

This strategy lead to a non-uniform loading dur- 
ing integration but leads to less message passing. 

Our experience has shown that Strategy II 
seems to work well on different computer platforms: 
both massively parallel computers as well as hetero- 
geneous cluster of workstations. So in the present ver- 
sion of the code, we have opted to implement Strategy 
II over Strategy I. 

16 DETAILS OF PARALLEL 
PERFORMANCE 

The details of the combined parallel perfor- 
mance of the CFD, EUPDF, and LSPRAY codes in- 
volving several different cases can be found in Refs. 


[10-15]. Here, we only summarize briefly the parallel- 
performance results for two different cases. One is 
a 3D test case and more details on this case can be 
found in the reference [13]. For this case, the calcu- 
lations were performed on a computational grid com- 
prising of 8430 tetrahedral elements and 100 Monte 
Carlo PDF particles per cell. The computations were 
performed on one of the NASA Ames Research Cen- 
ter’s parallel computer platforms called Turing which 
is a SGI Origin work-station with 24 PEs (Proces- 
sor Elements). Table 11 summarizes the CPU times 
per cycle taken by the EUPDF, LSPRAY, and CFD 
solvers vs the number of PEs. Both the CFD and 
PDF solvers show good parallel performance with an 
increase in the number of processors but for the spray 
solver it shows reasonable parallel performance. 

Next, we would like to summarize the results 
from [13, 38] showing only the results of spray compu- 
tations. The results are summarized in Table 12. The 
computations were performed on Turing at NASA 
ARC (Ames Research Center) - it is a SGI Origin 
work-station with a maximum of 24 PEs (Processor 
Elements). The computations made use of an un- 
structured grid with a mesh size of 3600 tetrahedral 
elements. And it made use of about 2,695 Lagrangian 
particles for the spray computations and one hundred 
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Monte Carlo particles per element for the PDF com- 
putations. In a given cycle of one global time step, 
dtgi, the spray equations were advanced over one hun- 
dred time steps as given by dt g i/dt m i = 100. The 
first row of Table 12 summarizes the CPU times per 
PE per cycle taken by the spray solver vs the num- 
ber of PEs. As expected, the CPU time goes down 
with an increase in the number of processors. In an 
ideal case, one would expect an inverse reduction in 
cpu time with an increase in the number of proces- 
sors. Here, we don’t get such an ideal performance 
because of the resulting non-uniform distribution of 
spray particles, between various participating proces- 
sors, from the implementation of Stratergy II. To get 
an idea of the spray particle distribution, we have tab- 
ulated the maximum and minimum number of parti- 
cles found between various processors. When we go 
from 1 to 2 PEs, 2097 particles are assigned to one 
and the rest to the second. With four they are dis- 
tributed between 1165 and 118, with eight between 
623 and 14, and with sixteen from 312 to 0. The 
results clearly show that the reduction in the CPU 
time varies almost linearly with the reduction in the 
number of maximum particles. 

17 SUMMARY OF SOME VALIDATION 

CASES INVOLVING BOTH REACTING 
AND NON-REACTING SPRAY 
COMPUTATIONS 

A total of the following five cases were validated: 

1. A reacting methanol spray with no-swirl. 

2. A non-reacting methanol spray with no-swirl. 

3. A confined swirl-stabilized n-heptane reacting 
spray. 

4. An unconfined swirl-stabilized n-heptane react- 
ing spray. 

5. A confined swirl-stabilized kerosene reacting 
spray. 

The experimental data for the first two cases 
was provided by McDonell & Samuelsen from the 
University of California at Irvine [38] . Both the cases 
are without swirl; one is a reacting case and the other 
is non-reacting. The data for the third and fourth 
cases was provided by Bulzan from the NASA Glenn 
Research Center [39-40]. Both the cases are swirl- 
stabilized reacting cases, one is an unconfined flame 


and the other is confined. The data for the last case 
was provided by El Banhawy & Whitelaw from Im- 
perial College [4]. It is a confined swirl-stabilized 
kerosene spray flame. The first three cases made 
use of unstructured grids and the last two structured 
grids. 

Here, we would like to provide a brief summary 
of the validation cases but a detailed presentation of 
the results and discussion can be found elsewhere in 
the papers [10-13]. The comparisons involved both 
gas and drop velocities, drop size distributions, drop 
spreading rates, and gas temperatures. The results 
were in reasonable agreement with the available ex- 
perimental data. The comparisons also involved the 
results obtained from the use of the Monte Carlo PDF 
method as well as those obtained from a conventional 
CFD solution without the Monte Carlo PDF method. 
For the first case of McDonell & Samuelsen’s reacting 
spray flame, the detailed comparisons clearly high- 
lighted the importance of chemistry /turbulence in- 
teractions in the modeling of reacting sprays [13]. 
The results from the PDF and non-PDF methods 
were found to be markedly different with the PDF 
solution providing a better approximation to the re- 
ported experimental data. The PDF solution showed 
that most of the combustion occured in a predom- 
inantly diffusion-type of flame environment and the 
rest occuring in a predominantly premixed-type of 
flame environment. However, the non-PDF predic- 
tions showed incorrectly that most of the combustion 
occured in a predominantly vaporization-controlled 
regime. The Monte Carlo temperature distribution 
showed that the functional form of the PDF for the 
temperature fluctuations varied substantially from 
point to point. The results brought to the fore some 
of the deficiencies associated with the use of assumed- 
shape PDF methods in spray computations. 

18 CONCLUDING REMARKS 

• This manual provides a complete description of 
LSPRAY-III - a Lagrangian spray solver devel- 
oped for application with parallel computing and 
unstructured grids. 

• It facilitates the calculation of the multi- 
component liquid sprays with variable properties 
valid over a wide range of low pressure condi- 
tions. 

• It provides the user with a basic understanding 
of the the spray formulation and the LSPRAY- 
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Ill code structure, and complete details on how 
to couple the spray code to any other flow code. 

• The basic structure adopted for the grid repre- 
sentation and parallelization for the gas side of 
the flow computations follows the guidelines es- 
tablished for NCC. 

• Also, we have extended the joint scalar Monte 
Carlo PDF method to two-phase flows and, 
thereby, demonstrating the importance of chem- 
istry/turbulence interactions in the modeling of 
reacting sprays. 


• Based on the validation studies involving several 
confined and unconfined spray flames, the results 
were found to be encouraging in terms of their 
ability to capture the overall structure of a spray 
flame. 

• The source code of LSPRAY-III will be available 
with NCC as a complete package. 
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APPENDIX 

DETAILS OF THE PRIMARY SPRAY 
ATOMIZATION MODELS 

1. Sheet Breakup Primary Atomization 
Model (Likely Application: Pressure Swirl 
Atomizer) 

Here, we summarize the details of the sheet 
breakup model taken from Schmidt et al [48]. The 
growth of an infinitesimal disturbance as given by 

p = r; 0 exp(ikx + cot) (1) 

was analyzed based on a linear stability analysis of a 
two-dimensional, viscous, incompressible liquid sheet 
of thickness 2 h which moves through an inviscid, in- 
compressible gas medium. This was analyzed in a co- 
ordinate system moving with the sheet with a relative 
velocity of U where p 0 is the initial wave amplitude, 
k (= 27t/A) is the wave number, and u> = to r + icoi 
is the complex growth rate. The most unstable dis- 
turbance responsible for the sheet breakup is denoted 
by 12. 

Based on the linearized liquid and gas continu- 
ity and momentum equations subject to the linearized 
boundary conditions at the gas and liquid interfaces, 
a sinuous mode dispersion relation was obtained by 
[49], 


uj 2 [tanh(kh) + Q] + cj [ 4 i 'ik 2 tanh(kh) + 2 iQkU]-\- 
4 vf k 4 tanh(kh) — 4 v 2 k 3 l tanh(lk ) — QU 2 k 2 + = 0 ( 2 ) 

where Q = p g /pi , l 2 = k 2 + uj/vi, and U is the 
relative velocity between the liquid and gas. Inviscid 
analysis also indicates that for low gas Weber num- 
ber (We = p g hU 2 /a) flows, long waves tend to grow 
leading to liquid sheet breakup but for higher Weber 
numbers, short waves produce a maximum growth 
rate followed by breakup. The critical Weber number 
that leads to the transition from the long wavelength 
regime to the short wavelength regime was shown to 
be We g = 27/16. For most modern fuel injection 
systems, the film Weber number is well above this 
critical limit. The growth rate for the sinuous mode, 
ay, based on an order of magnitude analysis of the 
dispersion relation yields, 

2uik 2 tanh{kK) 

ClJ'p 1 — 

tanh(kh) + Q 


4 vf 


k 4 tanh 2 (kh) — Q 2 U 2 k 2 — [tanh(kh) + Q]{—QU 2 k 2 + <rk 3 / pi) 
tanh(kh) + Q 


(3) 


For short waves in the limit of Q « 1 for high-speed 
sheets, it yields 


ay 


-2utk 2 + \Uvfk A + QU 2 k 2 


ak 3 

Pi 


( 4 ) 


Following Dombrowski & Johns [51], the sheet 
disintegration leads to the formation of ligaments 
once the unstable waves reach a critical amplitude 
and Eq. (4) shows that the growth rate of short 
waves is independent of the sheet thickness. The cor- 
responding breakup time r and the breakup length L 
are given by: 


p b = p 0 exp(Clr ) => r = ^.In— (5) 

12 p 0 

L = Vt = TjMj) (6) 

i 6 I Iq 

where 12 is the maximum growth rate as determined 
by Eq. (4), the term ln( — ) has an assigned value of 
12 as suggested by Dombrowski & Hooper [50], and 
V is the absolute velocity of the liquid. 

The initial diameter of the ligaments is derived 
from a mass balance relationship. For long waves, it 
is assumed that the ligaments are formed from tears 
in the sheet once per wavelength and the resulting 
diameter is given by, 


d-L = 



( 7 ) 


where K sb is the wave number corresponding to the 
maximum growth rate 12 as obtained from Eq. (3) 
and the film thickness, h , is calculated from the 
breakup length, L , the radial distance from the cen- 
terline to the mid-line of the liquid sheet at the at- 
omizer exit, r 0 , and the spray angle, 0, as follows: 
h = 2 T p,v(r 0 +L«n( 8 / 2 )) - For short waves, the lig- 
ament diameter is independent of the liquid sheet 
thickness and is assumed to be proportional to the 
wave length associated with the maximum growth 
rate 12 as follows: cIl = 2 j/ c / - where the ligament 
constant, Cl, is equal to 0.5. 

For both long and short waves, Dombrowski & 
Johns [51] developed a linear stability analysis for 
the further breakup from ligaments to droplets based 
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on the Weber’s analysis on capillary instability. The 
analysis shows that the breakup occurs when the 
amplitude of the unstable waves nears the radius of 
the ligament. And the corresponding most unstable 
wavenumber, Ku,, is given by: 


Ku,d,L 


1 + 1-1 

2 2 pi<rd L ) 


(8) 


This analysis thus yields the most probable droplet 
size based on a simple mass balance calculation where 
dl = 3Trd 2 L /K Lb . 


1.1 Application to pressure swirl atomization 


For the pressure swirl atomizer, the initial in- 
jector exit velocity and liquid sheet thickness are cal- 
culated following the approach of Schmidt et al [48] . 
It assumes a uniform velocity profile for the initial 
liquid velocity, V, as given by, 


V = 


max{ 0.7, 


4m 

7T d 2 picos6 



(9) 


where m and 6 are the measured mass flow rate and 
spray angle, respectively, d 0 is the injector exit di- 
ameter, and A p is the pressure drop in the injector. 
Once V is known, the corresponding axial component 
of the sheet velocity is calculated via u = V cosd. 
And the initial sheet thickness h Q is calculated from 
the conservation of mass: 


m = irpiuh 0 (d 0 - h 0 ) (10) 

At the point of primary breakup, the actual drop size 
is chosen from a Rosin- Rammler distribution with the 
mean size as given by d b of the sheet breakup model. 
Further movement of the droplets is tracked by mak- 
ing use of a Lagrangian formulation. 

2. Blob Jet Primary Atomization Model 
(Likely Application: Single-Orifice Nozzles) 

Here we summarize the details of the blob 
jet primary atomization model taken from Reitz & 
Bracco [41] and Reitz [52,45]. It applies for a cylin- 
drical liquid jet issuing into an incompressible gas. 
The following dispersion relation was obtained based 
on the stability analysis of a cylindrical liquid surface 
subjected to linear perturbations: 


2 o ,2 r I'i(k a ) 2 kl I\(ka) I[(la) 1 

W 1 ^ I 0 (ka) fc 2 + l 2 I 0 {ka ) I 0 (la ) 1 


ak 

Pia 2 


(1 -fcV){ 


l 2 — k 2 I\{ka ) 
l 2 + k 2 h 0 (ka) 




( 11 ) 


where Jo, I\, and Kq, K\ are the modified Bessel 
functions of the first and the second kinds. 

Reitz [52,45] generated curve-fits of numerical 
solutions to Eq. (11) for the maximum growth rate 
(cj = fl) and the corresponding wavelength (A = A): 


(1 + 0.45Z°' 5 )(1 + 0.4T 0 ' 7 ) , 

9 '° (1 + 0.87Wei-67)0.6 (12) 

0.34 + 0.38 Wef 5 

(13) 

(1 + Z)(l + 1.4T 0 - 6 ) ’ 

where Z = T = ZWe° g 5 , Wei = We g = 

pg ^ a , and Re i = A core region is predicted with 
the blob injection method because there is a region of 
large discrete liquid parcels near the nozzle. Based on 
the jet stability theory, new drops are formed from a 
parent drop or blob. It is assumed that small droplets 
(with radius, r) are formed from the parent drops 
(with radius, a) with drop size proportional to the 
wavelength of the fastest-growing or most-unstable 
wave, 


n {^-} 0 - 5 


r = B 0 A (if B a A < a) or 

r = min{(3Tra 2 U /2fl) 0 ' 33 , (3a 2 A/4) 0 ' 33 } 

(if B 0 A > a, one time only) (14) 

where B a = 0.61 according to Reitz [52,45]. In the 
above, it is assumed for the (B a A < a) condition 
that small droplets are formed with the dropsize pro- 
portional to the wavelength of the fastest growing 
mode and the second (B 0 A > a) condition applies to 
drops larger than the jet and it assumes that the jet 
disturbance has frequency 12/27 r ( a drop is formed 
each wave period) or that the drop size is determined 
from the volume of liquid contained under one sur- 
face wave. And the rate of change of droplet radius 
due to breakup is given by 


da/dt = —(a — r)/r (r < a) (15) 

where r is the breakup time, r = 3.726-Bia/Afl, and 
the value for the breakup time constant, B i, depends 
on the injector characteristics and its value ranges 
between 1.732 to 40. And aft = t 0 ) = a 0 is the 
initial drop radius at time, t Q . 
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After the breakup, a new parcel containing 
product drops of size, r, is created and added to the 
computations. This was done as long as the mass of 
the liquid removed from the parent (/9/47r(a 3 — a 3 )/3) 
reached or exceeded 3% of the average injected parcel 
mass and if the number of product drops is greater 
than or equal to the number of parent drops [45]. 
While waiting for sufficient product drops to accu- 
mulate, the parent drop number was adjusted so that 
Na 3 = _/V 0 a 3 but the parent drop number, N a , was 
then restored following the creation of the new prod- 
uct parcel [45]. 

In the case of (B a A > a), the parent parcel 
was replaced by a new parcel containing drops with 
size given by Eq. 14 after a time equal to r (with 
N = N Q a 3 /r 3 ) [45]. This breakup procedure was 
allowed only once for each injected parcel [45]. 

Validation of the model for a single hole orifice 
in a typical diesel engine was demonstrated by Reitz 
and Diwaker [46] and Reitz [52,45]. 


Here we summarize the details of the BLS at- 
omization model taken from Khosla and Crocker [47] . 
In this model both surface shear breakup and col- 
umn breakup modes are included. Before column 
breakup, fragments may be formed due to boundary 
layer stripping depending on the local Weber num- 
ber and q (= piuf / p g v? g ). When the jet reaches the 
column breakup time, the entire jet breaks into frag- 
ments. It also allows for further breakup of the frag- 
ments based on a modified boundary layer stripping. 
And it is followed by a final breakup step based on the 
Rayleigh-Taylor secondary droplet breakup model. 

Fragments are stripped from the liquid column 
if the We g satisfies the following criteria: 

We g > 50Re g /2 9 " 1/0 - 81 and 

We g > 15 (18) 

where 


3. Air Blast Primary Atomization Model 
(Likely Application: Air Blast Atomizers) 

The air blast atomization model is essentially 
based on the idea of pressure-swirl atomization model 
(Section 1.1) as the primary atomization of an air 
blast injector is based on the aerodynamic analysis 
involving the Kelvin-Helmholtz instability of a liq- 
uid jet in an incompressible gas. It, however, differs 
from the pressure-swirl atomization model in the de- 
termination of the initial sheet velocity and thickness, 
which are given by 


V sheet = aV x + (l-a)V g (16) 


S 


r [l 



rh 


7T r 2 pV sheet 


(17) 


where a has a value of 0.12 to 1 depending on the 
fuel filmer characteristics, r is the radius of the fuel 
filmer, and rh is the fuel flow rate. The continuous 
annular sheet is represented by a finite number of 
point injectors located randomly along the circular 
ring of the liquid sheet. 

4. Modified BLS Primary Atomization 
Model (Likely Application: Liquid Jet in a 
Cross Flow) 


Re g — PgdjUg l p g 

we g = u 2 g djp g /ai 

where Re g and We g are based on the gas velocity 
component normal to the liquid jet direction instead 
of the relative velocity. If column stripping does oc- 
cur, the amount of mass removed from the column is 
given by 


where, 


h 


7 rd 


M shed = -Trdpi—UreiAa\/ — At 


A = 


a = 


r Pjl~\ V 3 rAhn V 3 
PZ MZ 

8 Pl 1 1/2 


r_UW_n 

L 3 AUrelpR 


t* = 


Hf. 


(19) 


( 20 ) 

(21) 


(22) 


where is the liquid column drop lifetime. The ad- 
dition of the factor t*/tb causes the shedding rate to 
increase essentially linearly with distance away from 
the injection location. This accounts for the lack of 
shedding close to the injection location and the sub- 
sequent buildup of shedding over the life of the liquid 
column. The shed drop SMD Sauter Mean Diameter) 
is given by 
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SMD = 3A-d 1/2 [— ] 1/4 [_^_] 1/2 (23) 

t* Pa Ur el Pi 

The t b /t* factor is included with the effect of 
producing smaller drops near the injection location. 
However, t b /t* is limited to a minimum value of 2.5 
which is never exceeded for some cases. The amount 
of mass shed is tracked and 10 new parcels are cre- 
ated when the cumulative shed mass after a time step 
exceeds 1% of the mass of the parcel. Each parcel is 
allocated an equal amount of the shed mass and the 
size for each new parcel is selected randomly from a 
uniform distribution between 0.4 and 1.6 times the 
mass mean diameter, MMD. And the drop velocities 
were given by 

u d = u P + 0.3(RND — 0.2) (u g — u p ) (24) 


Vd = v p + +0.3(RND — 0.2)(v g - v p ) (25) 


w d = w p + 0.25 (RND — 0.5) (u re i — w p ) (26) 


the same as Eqs. (24) and (25). The lateral velocity 
component is given by 

Wd = w p + 0A(RND — 0.5 )(u re i — w p ) (28) 

Fragments are further broken into small drops 
according to a modified version of the boundary layer 
stripping model based on the following criteria, 

We g > V Re and 

We g > 15 (29) 

where 

we g = u 2 rel d d p g /<Ji 

Re g = PgddU re i/ p g 

The criteria are generally the same as for column 
stripping except that the dependence on q is not 
needed. Also note that We and Re are now deter- 
mined using the relative velocity instead of the cross 
flow velocity. Again, the mass shed from a fragment 
in a time step and the SMD are given by 


where RND is a random number. Column stripping 
occurs, assuming the above criteria are met, until the 
column breakup time is exceeded. Parcels created 
through the column stripping mechanism are con- 
sidered to be fragments which may undergo further 
breakup as discussed below. First, though, the col- 
umn breakup mechanism, which also produces frag- 
ments, is described. The liquid jet column breakup 
time is given by 

t b = A b We 0 62 t* (27) 

The constant, A b , has a value of 25. Since the present 
model includes a fragment breakup step after the col- 
umn breakup, the We dependence for the breakup 
onset was retained. 

After the column breakup time is reached, the 
column is broken into 18 new parcels with MMD = 
0.45dj. The new parcels are also designated as frag- 
ments. The size of the fragments still tend to be 
large, so the ultimate drop size from the primary 
breakup process is mostly determined from the frag- 
ment breakup process. The size distribution is the 
same as described above for the column stripping. 
The cross flow and normal velocity components are 


M shed = 1.2ndpiu re iAay —A t (30) 

SMD = 3.6d 1/2 [—] 1/4 [— ^ — 1 1/2 (31) 

where A and a are given by Eqs. (20) and (21), re- 
spectively. The new droplet velocities are given by 
Eqs. (24), (25), and (28). The broken fragments pro- 
duce 3 new parcels with size distribution the same as 
described above when the shed mass from the frag- 
ment exceeds 20% of the fragment mass. A fragment 
can continue to breakup until it no longer meets the 
criteria of Eq. (29) or until its size is lower than the 
newly created drops. Once the fragment breakup pro- 
cess is complete, drops may breakup further based on 
a Rayleigh-Taylor secondary breakup method. 

Khosla and Crocker [47] applied the model to 
predict the properties of Jet A-l kerosene fuel in- 
jected into a cross-flowing air stream. 

DETAILS OF THE SECONDARY 
DROPLET BREAKUP MODELS 
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Recent studies by Reitz et al [55, 52] have ex- 
amined the breakup of single droplets moving in 
a transverse, high-velocity air jet. The high-speed 
photography provided new insights into the details 
of the breakup mechanism of a single drop. The 
droplet breakup regimes are classified as bag, strip- 
ping (shear) and catastrophic (surface wave) based 
on an increasing size of Weber number. In the bag 
breakup mode (at low Weber number), the drop is 
flattened by the aerodynamic pressure, then turned 
inside out, forming a thin hollow bag which is tied to- 
gether with a circular belt-like structure on the wind- 
ward side. The bag eventually bursts into smaller 
liquid fragments, whereas the belt decays into larger 
ligaments and droplets. In the stripping regime thin 
sheets or ligaments of fluid are continuously shed from 
the periphery of the distorting parent drop as a con- 
sequence of a K-H instability, causing these sheets to 
disintegrate into tiny droplets. This process always 
leaves a coherent residual parent drop. The catas- 
trophic drop breakup takes place in two stages lead- 
ing to a collection of larger and tiny product droplets: 
Large amplitude long-wavelength waves caused by 
drop deceleration induce a R-T instability on the 
flattened drop which leads to a breakup into large 
product droplets, while at the same time short sur- 
face waves induce a K-H instability on the windward 
side of the parent drop resulting in a collection of 
much smaller product droplets. In diesel or other 
high-pressure gas-turbine sprays the droplets span a 
wide range of velocities and hence Weber numbers, 
and thus it is expected that all three droplet breakup 
mechanisms are relevant in the breakup modeling. 

In what follows, we provide some details of 
the secondary droplet breakup models contained in 
the CFDRC/UW atomization module: (1) Rayleigh- 
Taylor, (2) TAB, and (3) ETAB. These details are 
taken from [52-56, 42-44]. 

1. Rayleigh- Taylor Secondary Droplet 
Breakup Model 

Here we summarize the Rayleigh-Taylor sec- 
ondary breakup model developed by Patterson and 
Reitz [57]. It is based on the analysis developed 
by Taylor [53,57] that accounts for the disturbances 
caused by droplet deceleration. In the Rayleigh- 
Taylor breakup mechanism, the breakup wavelength, 
A, is given by 

A = Q.2-K^Za/\ud\{pi - p g ) (32) 


where ii,j is the drop deceleration (= | Cd ^ a J 7 , Cd is 
the aerodynamic drag coefficient, U is the drop rela- 
tive velocity, & d is the drop diameter) . Furthermore, 
the breakup wavelength A is limited by 


A = max(0.8d, A)[l + 0.2(RND — 0.5)] (33) 

and the breakup time, tb,RTi is given by 


tb,RT = Cfreq > 0.5^(pi + Pg)( 


\Ud\(pi~Pg) 


1.5 


(34) 

However, the value assigned for the constant, Cf req , 
depends on the droplet classification - parent, prod- 
uct, or default. For more details, one can refer to 
Patterson and Reitz [21]. After the breakup, no new 
drop parcels are created and there is no change in ve- 
locities between the parent and product drops. How- 
ever, the drop number in a given parcel changes to 
n product as given by n parent (d parent /A) 3 due to the 
change in the sizes between the parent and product 
drops. 


2. The TAB Secondary Droplet Breakup 
Model (Likely Application: 

Lenticular-Shaped Droplet Deformations) 

In an attempt to provide a description of the 
droplet and jet disintegration in the modeling of 
diesel sprays, O’Rourke & Amsden introduced the 
TAB model [42] . Here we summarize the TAB model 
taken from [43-44]. 

The TAB breakup model is based on Taylor’s 
analogy between an oscillating, distorting drop and a 
spring-mass system. A detailed analysis of this model 
together with a discussion of its numerical implemen- 
tation can be found in [42 & 58]. In this model, the 
drop motion is governed by a linear ordinary differ- 
ential equation for a forced, damped harmonic oscil- 
lator. The forcing term is given by the aerodynamic 
droplet-gas interaction, the damping is due to the 
liquid viscosity and the restoring force is supplied by 
the surface tension. The parameters and constants 
have been determined partly from theoretical consid- 
erations and partly from experimental observations. 

The TAB model describes the distortion of the 
drop by the deformation parameter, y = 2x/a, 

where x denotes the increase in the radius increase 
of the equator from its equilibrium position and a 
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is the drop radius. The equation for the distortion 
parameter y is given by 


V + 


5 m 

Pia 2 


y + 



2p g \u\ 2 

3 pia 2 


(35) 


Assuming a constant relative drop-gas velocity, 
U (which is satisfied in the numerical solution process 
during a given time step), the solution to Eq. (35) is 
given by 


. , We 

M = 12 + 


=i n We-, 

e *d ( [y( 0 ) — J cosojt 


-y 

L u; 


2/(0) - if - 

utd 


smio 


l ) 


(36) 


Here we provide some details of the ETAB 
model taken from Tanner [43-44]. The ETAB model 
is based on the following modifications to the stan- 
dard TAB model: (1) the droplet disintegration is 
modeled via an exponential law which relates the 
mean product droplet size to the breakup time of the 
parent drop; and (2) an energy balance consideration 
between the parent and product droplets yields an 
expression for the normal velocity component of the 
product droplet. 

When the breakup condition of We > 
W e cr it = 6 and y(t) = 1 is met, then the parent 
drop breaks up into a collection of product droplets, 
subject to a size distribution function which, in gen- 
eral, depends on the breakup mechanism. In the 
ETAB model, the rate of product droplet generation, 
dn(t)/dt, is given by 


where Ife = p g aU 2 /a and 


td 


2pia 2 
5 pi 


(37) 


8a 1 

Pia 3 t 2 


(38) 


In this model, it is assumed that a necessary 
condition for drop breakup is reached when We > 
We cr it . And the value for the critical Weber num- 
ber is determined experimentally to be 6. For an 
inviscid liquid with initial deformation conditions 
2/(0) = 2/(0) = 0, the solution to Eq. (35) 

reduces to y(t) = We( 1 — cosw 0 t)/ 12, where 

u) 2 = [44]. Consequently, the breakup occurs 

when y(t) > 1. The drop size after breakup is de- 
termined by an energy balance equation between the 
parent and the product droplets which equates the 
surface energies of parent drops with the energies of 
product drops due to oscillation and distortion. Also, 
the product droplets are initially equipped with the 
additional deformation velocity x = ay/2, which 

acts normal to the path of the parent drop and is 
responsible for the formation of the spray angle. 

One major advantage of this model is that it is 
based on a simple linear equation and it can be used 
effectively to describe the lenticular-shaped droplet 
deformations as observed in the experiments of [56 & 
54], 


3. The ETAB Secondary Droplet Breakup 
Model (Likely Application: High-Pressure 
Diesel Engine) 




(39) 


where n{t) = mo/rh(t ) and mo is the mass of the 
parent drop and to the mean mass of the product 
droplet distribution. Utilizing the fact that dn/dt = 
— ( mo /in 2 ) (dm /dt ), leads to the breakup law which 
relates the product drop size to the breakup time as 
determined by the TAB model. 


— = -3K b m (40) 

The breakup constant K b depends on the 
breakup regime and is given by parent drop prop- 
erties only. Bag breakup occurs if We = Wet and 
stripping breakup happens if We > We t . And it is 
given by 


I\ b = k\LO if We < We t or 

K b = k 2 uVWe if We > We t (41) 

The values for Wet, fci and k 2 have been determined 
experimentally and has been set to k\ « fc 2 = 1/4.5 
and We t = 80. 

In this model, a uniform product droplet size 
distribution is assumed. It is also noted that the 
choice of uniform distribution is not expected to be 
realistic but may produce good approximations when 
averaged over many drop breakups, because parent 
drops of different sizes and Weber numbers will in 
general yield a wide range of duct droplet sizes. With 
this assumption, Eq. (40) becomes 

- = e~ Kbt (42) 

a 
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where a and r are the radii of the parent and product 
drops, respectively. 

After breakup of a parent drop, the initial de- 
formation parameters of the product droplets are set 
to y( 0) = y(0) = 0. Also, the product droplets are 
initially supplied with a velocity component perpen- 
dicular to the path of the parent drop with a value 
vt = Ax, where A is a constant determined from the 
following energy balance consideration. The energy 
of the parent drop is the sum of the surface tension en- 
ergy and the droplet deformation energy. The second 
is computed as the product of the aerodynamic drag 
and the drop deformation at the stagnation point, 
estimated to be 5a/9. This leads to 

Eparent = 4 ttcto 2 + 57rc d p g a 3 1 U\ 2 /18 (43) 

And the energy of the product droplets in the frame 
of reference of the parent drop is given by 


A 2 = 3[1 - a/rsMR + 5c d W e / 72}co 2 /y 2 (45) 

where w 2 = 

Pia 6 

Tanner [43-44] analysis yields an approximate 
value of 0.69 for A showing that only 70% of the par- 
ent drop deformation velocity goes into the normal 
velocity component of the product droplets, where as 
it is 100% in the standard TAB model. 

Also, the characteristic time, r (= ^ ), for 
breakup in Eq. (42) for an inviscid liquid (/q = 

0) is given by a\ J if We < We t or 
a2 y^PgW\ if We > We t , where the suggested val- 
ues are for a\ = (-\/8A;i) -1 and c*2 = (\J 8^2) _ 1 - 

The application of the ETAB model in the sim- 
ulation of a high pressure liquid jet breakup can be 
found in [43-44]. 


Eproduct = 47 Taa 3 /r S MR + A 2 npia 5 y 2 /6 (44) 

where the Sauter mean radius, tsmr , enters via the 
relation r 2 = o 3 /tsmr- From Eqs. (43) and (44) 
one obtains the relation 


NASA/CR— 2008-215290 


45 



REFERENCES 

1. W.A. Sirignano, “Fluid Dynamics of Sprays,” 
Journal of Fluids Engineering, vol. 115, no. 3, 
pp. 345-378, September 1993. 

2. Crowe, C.T., “Numerical Models for Dilute Gas- 
Particle Flows,” Journal of Fluids, Vol. 104, pp. 
297-301, 1982. 

3. C.T. Crowe, M.P. Sharma, and D.E. Stock, The 
Particle-Source-in Cell (PSI-CELL) Model for 
Gas-Droplet Flows, J. Fluids Eng., vol. 99, pp. 
325, 1977. 

4. Y. El Banhawy and J.H. Whitelaw, Calculation 
of the Flow Properties of a Confined Kerosene- 
Spray Flame, AIAA J., vol. 18, no. 12, pp. 
1503-1510, 1980. 

5. Faeth, G.M., “Mixing, transport, and Combus- 
tion in Sprays,” Progress Energy Combustion 
Science, Vol. 13, pp. 293-345, 1987. 

6. Raju, M.S., and Sirignano, W.A., “Spray Com- 
putations in a Centerbody Combustor,” Pro- 
ceedings of the 1987 ASME-JSME Thermal En- 
gineering Joint Conference, Vol. 1, pp. 61-71, 
Honolulu, HI, March 1987. Also see Journal of 
Engineering for Gas Turbines and Power, Vol. 1, 
No. 4, pp. 710-718, October 1989. 

7. Raju M.S., and Sirignano, W.A., “Multi- 
Component Spray Computations in a Modified 
Centerbody Combustor,” Journal of Propulsion 
and Power, Vol. 6, No. 2, March-April 1990. 

8. Raju, M.S., “AGNI-3D: A Computer Code for 
the Three-Dimensional Modeling of a Wankel 
Engine,” Computers in Engine Technology, Pro- 
ceedings IMechE, London, United Kingdom, pp. 
27-37, 1991. 

9. Raju, M.S., “Heat Transfer and Performance 
Characteristics of a Dual-Ignition Wankel En- 
gine,” Journal of Engines, the 1992 SAE Trans- 
actions, Section 3, pp. 466-509. 

10. M.S. Raju, Application of Scalar Monte Carlo 
Probability Density Function Method For Tur- 
bulent Spray Flames, Numerical Heat Transfer, 
Part A, vol.' 30, pp. 753-777, 1996. 


11. M.S. Raju, Current Status of the Use of Par- 
allel Computing in Turbulent Reacting Flows: 
Computations Involving Sprays, Scalar Monte 
Carlo Probability Density Function & Unstruc- 
tured Grids, Advances in Numerical Heat Trans- 
fer, vol. 2, ch. 8, pp. 259-287, 2000. 

12. M.S. Raju, Scalar Monte Carlo PDF Com- 
putations of Spray Flames on Unstructured 
Grids With Parallel Computing, Numerical Heat 
Transfer, Part B, No. 2, Vol. 35, pp. 185-209, 
March 1999. 

13. M.S. Raju, On the Importance of Chem- 
istry/Turbulence Interactions in Spray Compu- 
tations, Numerical Heat Transfer, Part B: Fun- 
damentals, No. 5, Vol. 41, pp. 409-432, 2002. 

14. M.S. Raju, LSPRAY - A Lagrangian Spray 
Solver - User’s Manual, NASA/CR.-97-206240, 
NASA Lewis Research Center, Cleveland, Ohio, 
November 1997. 

15. M.S. Raju, EUPDF - An Eulerian-Based Monte 
Carlo Probability Density Function (PDF) 
Solver - User’s Manual, NASA/CR.-1998-20401, 
NASA Lewis Research Center, Cleveland, Ohio, 
April, 1998. 

16. Reid, R.C., Prausnitz, J.M., and Poling, B.E., 
“The Properties of Gases and Liquids,” 4th Edi- 
tion, McGraw-Hill Book Company, 1987. 

17. Prausnitz, J.M., and Chueh, P.L., “Computer 
Calculation For High pressure Vapor-Liquid 
Equilibria,” Prentice-Hall, Englewood Cliffs, 
N.J., 1968. 

18. Aggarwal, S.K., Shu, Z., Mongia, H., and Hura, 
H., “Multicomponent Fuel Effects on the Vapor- 
ization of a Surrogate Single-Component Fuel 
Droplet,” AIAA Paper 98-0157, 36th Aerospace 
Sciences Meeting, Reno, Nevada, Jan. 12-15, 
1998. 

19. Aggarwal, S.K., Shu, Z., Mongia, H., and Hura, 
H., “Multicomponent Single-Component Fuel 
Droplet Evaporation Under High Pressure Con- 
ditions,” 34th AIAA/ASME/SAE/ASEE Joint 
Propulsion Conference & Exhibit, Cleveland, 
OH, July 13-15, 1998. 

20. Bellan, J., “Supercritical (and Subcritical) fluid 
Behavior and Modelling: Drops, Streams, Shear 


NASA/CR— 2008-215290 


46 



and Mixing Layers, Jets, and sprays,” Prog. En- 
ergy Combust. Sci., Vol. 26, pp. 329-366, 2000. 

21. R. Ryder, CORSAIR User’s Manual: Version 
1.0, SID: Y965, Pratt and Whitney Engineer- 
ing, United Technologies Corporation, 25 Jan- 
uary 1993. 

22. N.S. Liu and R.M. Stubbs, Preview of 
National Combustion Code, AIAA 97-3114, 
AIAA/ASME/SAE/ASEE 33rd Joint Propul- 
sion Conference, Seattle, Wash., July 6-9, 1997. 

23. K.-H. Chen, A.T. Norris, A. Quealy, and 
N.-S. Liu, Benchmark Test Cases For the 
National Combustion Code, AIAA 98-3855, 
AIAA/ASME/SAE/ASEE 34th Joint Propul- 
sion Conference, Cleveland, Ohio, July 13-15, 
1998. 

24. McBride, S. Gordan, and M. Reno, “Coefficients 
for Calculating Thermodynamic and Transport 
Properties of Individual species,” NASA TM- 
4513, NASA Lewis Research Center, 1993. 

25. C.R. Wilke, “A Viscosity Equation for Gas Mix- 
ture,” Chem. Physics, Vol. 18, no. 4, pp. 517- 
519, April 1950. 

26. F.M. White, Viscous Flows, McGraw-Hill Inc., 

1974. 

27. K. Kuo, Principles of Combustion, Wiley and 
Sons, 1976. 

28. Peng, D.Y., and Robinson, D.B., Ind. Eng. 
Chem. Fundam., Vol. 15, pp. 59, 1976. 

29. Reichenberg, D., “The Viscosities of Gases at 
High Pressures,” Natl. Eng. Lab., Rept. Chem. 
38, East Kilbride, Glascow, Scotland, August 

1975. 

30. Stiel, L.I., and Thodos, G., AICHE J., Vol. 10, 
pp. 26, 1964. 

31. Takahashi, S., J. Chem. Eng. Japan, Vol. 7, pp. 
417, 1974. 

32. Yuen, M.C., and Chen, L.W., “On Drag of Evap- 
orating Droplets,” Combust. Sci. Technol., Vol. 
14, pp. 147-154, 1976. 

33. A. A. Amsden, P.J. O’Rourke, and T.D. Butler, 
“KIVA-II: A Computer program For Chemically 
Reactive Flows With Sprays,” LA-11560-MS, 


Los Alamos National Laboratory, Los Alamos, 
New Mexico 87545. 

34. A.Y. Tong and W.A. Sirignano, Multi-compon- 
ent Transient Droplet Vaporization With Inter- 
nal Circulation: Integral Formulation and Ap- 
proximate Solution, Numerical Heat Transfer, 
vol. 10, pp. 253-278, 1986. 

35. R. Clift, J.R. Grace, and M.E. Weber, Bub- 
bles, Drops, and Particles, Academic, New York, 
1978. 

36. H. Schlichting, Boundary-Layer Theory, Mc- 
Graw-Hill Series in Mechanical Engineering: 
McGraw-Hill, Inc., New York, 1968. 

37. J.S. Ryan and S.K. Weeratunga, Parallel Com- 
putation of 3-D Navier-Stokes Flowfields for Su- 
personic Vehicles, AIAA 93-0064, AIAA 31st 
Aerospace Sciences Meeting, Reno, Nevada, 
1993. 

38. V.G. McDonell and G.S. Samuelson, An Ex- 
perimental Data Base for the Computational 
Fluid Dynamics of Reacting and Nonreacting 
Methanol Sprays, J. Fluids Engineering, vol. 
117, pp. 145-153, 1995. 

39. D.L. Bulzan, Structure of a Swirl-Stabilized, 
Combusting Spray, NASA Technical Memoran- 
dum: NASA TM-106724, Lewis Research Cen- 
ter, Cleveland, Ohio, 1994. 

40. D.L. Bulzan, Velocity and Drop Size Measure- 
ments in a Confined, Swirl-Stabilized Combust- 
ing Spray, AIAA 96-3164, 32rd AIAA/ ASME/ 
SAE/ ASEE Joint Propulsion Conference, July 
01-03, 1996/Buena Vista, FL. 

41. Reitz, R.D., & Bracco, F.V., “Mechanism of At- 
omization of a Liquid Jet,” Phy. Fluids, Vol. 25, 
ch. 10, Oct., 1982. 

42. O’Rourke, P.J., and Amsden, A. A., “The TAB 
Method for Numerical Calculation of Spray 
Droplet Breakup,” SAE Technical Paper 872089, 
1987. 

43. Tanner, F.X., “Liquid Jet Atomization and 
Droplet Breakup Modeling of Non-evaporating 
Diesel Fuel sprays,” SAE Technical Paper 
970050, 1998. Also SAE 1997 Transactions: 
Journal of Engines, Vol. 106, Sec. 3, pp. 127- 
MO, 1998. 


NASA/CR— 2008-215290 


47 



44. Tanner, F.X., “A Cascade Atomization and 
Drop Breakup Model for the Simulation of High- 
Pressure Liquid Jets,” SAE Paper 2003-01-1044, 
2003. 

45. Reitz, R.D., “Modeling Atomization Processes 
in High-Pressure vaporizing Sprays,” Atomiza- 
tion and Spray Technology, Vol. 3, pp. 309-337, 
1987. 

46. Reitz, R.D., and Diwakar, R. “Structure of High- 
Pressure Fuel Sprays,” SAE paper 870598, 1987. 

47. Khosla, S., and Crocker, D.S., “CFD Model- 
ing of the Atomization of Plain Jets in Cross 
Flow for Gas Turbine Applications” ’ IGTI Turbo 
Expo: Combustion & Fuels, GT2004-54269, Vi- 
enna, Austria, June 2004. 

48. Schmidt, D.P., Nouar, I., Senecal, P.K., Hoff- 
man, J., Rutland, C.J., Martin. J., & Reitz, 
R.D., “Pressure-Swirl Atomization in the Near 
Field,” 1999 SAE Congress, SAE 1999-01-0496. 

49. Senecal, P.K., Schmidt, D.P., Nouar, I., Rutland, 
C.J., & Reitz, R.D., “Modeling High Speed Vis- 
cous Liquid Sheet Approximation,” 1998. 

50. Dombrowski, N., & Hooper, P.C., “The Effect 
of Ambient Density on Droplet Formation in 
Sprays,” Chem. Eng. Sci. vol. 17, pp. 291, 
1962'. 

51. Dombrowski, N., & Johns, W.R., “The Aerody- 
namic Instability and Disintegration of Liquid 
Sheets,” Chem. Eng. Sci. vol. 1718, pp. 203, 
1963. 

52. Reitz, R.D., “Computer Modeling of Sprays,” 
Spray Technology Short Course, Pittsburg, PA 
May 7, 1997. 

53. Taylor, G.I., “The instability of Surfaces When 
Accelerated in a Direction Perpendicular to 
Their Planes,” Proc. Royal Soc., A, Vol. 201, 
pp. 192-196, 1950. 

54. Hwang, S.S., Liu, Z., and Reitz, R.D., “Breakup 
Mechanisms and Drag Coefficients of High- 
Speed Vaporizing Liquid Drops,” Atomization 
and Sprays, Vol. 6, pp. 353-376, 1996. 

55. Liu, A.B., Mather, D., and Reitz, R.D., “Mod- 
eling the Effects of Droplet Drag and Breakup 
on Fuel sprays,” SAE Technical paper 930072, 
1993. 


56. Liu, A.B., and Reitz, R.D., “Mechanisms of Air- 
Assisted Liquid Atomization,” Atomization and 
Sprays, Vol. 3, pp. 55-75, 1993. 

57. Patterson, M.A., and Reitz, R.D., “Modeling the 
Effects of Fuel Characteristics on Diesel Engine 
Combustion and Emission,” SAE 980131,1998. 

58. Amsden, A. A., O’Rouke, P.J., and Butler, T.D., 
“KIVA II: A Computer Program for Chemi- 
cally Reactive Flows With Sprays,” Technical re- 
port LA-11560-MS, Los Alamos laboratory, May 
1989. 

59. Raju, M.S., “Our Computational Experience 
With the Integration of a Spray Atomization 
Module,” ICCES0520050720481, International 
Conference on Computational & Experimental 
Engineering Sciences, Chennai, India, Dec. 1- 
10, 2005. 

60. Raju, M.S., “Numerical Investigation of Various 
Atomization Models in the Modeling of a Spray 
Flame,” AIAA Paper 2006-0176, 44th Aerospace 
Sciences Meeting and Exhibit, Reno, Nevada, 
Jan. 9-12, 2006. 

61. Lucas, K.D., Tseng, C.C., Pourpoint, T.L., 
Lucht, R.P., and Anderson, W.E., “Imaging 
Flashing Injection of Acetone at Jet Engine Aug- 
mentor Conditions,” AIAA Paper 2007-1182, 
45th AIAA Aerospace Sciences Meeting and Ex- 
hibit, Reno, Nevada, Jan. 8-11, 2007. 

62. Zuo, B., Gomes, A.M., and Rutland, C.J., 
“Modeling Superheated Fuel Sprays and Vapor- 
ization,” Inj. J. Engine Research, vol. 1, no. 4, 
pp. 321-336. 

63. Schmehl, R, and Steelant, J., 
“Flash-Evaporation of Oxidizer During Start-Up 
of an Upper-Stage Rocket Engine,” AIAA Paper 
2003-5075, 39th AIAA/ASME SAE/ASEE Joint 
Propulsion Conference and Exhibit, Huntsville, 
Alabama, July 20-23 2003. 

64. Schmehl, R, and Steelant, J., “Evaluation of Oxi- 
dizer Temperature Drop in a Combustion Cham- 
ber,” 4th International Conference on Launcher 
Technology ” Space Launcher Liquid Propulsion, 
Liege, Belgium, December 3-6 2002. 

65. Aclachi, M., McDonnel, V.G., Tanaka, D., 
Senda, J, and Fujimoto, H., “Characteristics of 


NASA/CR— 2008-215290 


48 



Fuel Vapor Concentration Inside a Flash Boiling 
Spray,” SAE Paper 970871, 1997. 

66. Lefebvre, A., “Atomization and Sprays,” Hemi- 
sphere Publishing Company, New York, pp. 165- 
222, 1989. 

67. VanDerWege, B.A., Lounberry, T.H., and 
Hochgreb, S., “Numerical Modeling of Fuel 
Sprays in DISI engines Under Early-Injection 
Operating Conditions,” SAE Paper 2000-01- 
0273, 2000. 

68. Reitz, R.D., “A Photographic Study of Flash- 
Boiling Atomization,” Aerosol Science Technol- 
ogy, Vol, 12, pp. 561-569, 1990. 

69. R. Borghi, Turbulent Combustion Modeling, 
Prog. Energy Combust. Sci., vol. 14, pp. 245- 
292, 1988. 

70. S.B. Pope, PDF Methods for Turbulent Reactive 
Flows, Prg. Energy Combust. Sci., vol. 11, pp. 
119-192, 1985. 


71. S.M. Correa, Development and Assessment 
of Turbulence-Chemistry Models in Highly 
Strained Non-Premixed Flames, AFOSR/NA 
Contractor Report, 110 Duncan Avenue, Bolling 
AFB, DC 20332-0001, 31 October 1994. 

72. Yildiz, D., Rambaud, P., Van Beeck, J., Buch- 
lin, J.-M., “Characterization of Superheated Liq- 
uid Jet Atomization Phase Doppler Anemome- 
ter (PDA) and High-Speed Imaging,” Proceed- 
ings of FEDSM2006: 2006 ASME Joint U.S.- 
European Fluids Engineering Summer Meeting, 
July 17-20, 2006, Miami, Florida. 

73. Yildiz, D., Rambaud, P., Van Beeck, J., Buchlin, 
J.-M., “Evolution of the Spray Characteristics in 
Superheated Liquid Jet Atomization in Function 
of Initial Flow Conditions,” ICLASS-2006, Paper 
ID ICLASS06-122, Aug. 27-Sept 1, 2006, Kyoto, 
Japan. 


NASA/CR— 2008-215290 


49 



REPORT DOCUMENTATION PAGE 

Form Approved 
OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 
data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 
burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 
Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 
control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1 . REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 

01-07-2008 Final Contractor Report 

3. DATES COVERED (From - To) 

4. TITLE AND SUBTITLE 

LSPRAY-II1: A Lagrangian Spray Module 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

NNC06BA07B 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

Raju, M.S. 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

22 

5f. WORK UNIT NUMBER 

WBS 984754.02.07.03.19.02 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

ASRC Aerospace Corporation 
Glenn Research Center 
Cleveland, Ohio 44135 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

E- 165 64 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

10. SPONSORING/MONITORS 
ACRONYM(S) 

NASA 

11. SPONSORING/MONITORING 
REPORT NUMBER 

NASA/CR-2008-215290 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified-Unlimited 

Subject Categories: 07, 28, 34, and 61 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 301-621-0390 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

LSPRAY -III is a Lagrangian spray solver developed for application with parallel computing and unstructured grids. It is designed to be 
massively parallel and could easily be coupled with any existing gas-phase flow and/or Monte Carlo Probability Density Function (PDF) 
solvers. The solver accommodates the use of an unstructured mesh with mixed elements of either triangular, quadrilateral, and/or tetrahedral 
type for the gas flow grid representation. It is mainly designed to predict the flow, thermal and transport properties of a rapidly vaporizing 
spray because of its importance in aerospace application. The manual provides the user with an understanding of various models involved in 
the spray formulation, its code structure and solution algorithm, and various other issues related to parallelization and its coupling with other 
solvers. With the development of LSPRAY -III, we have advanced the state-of-the-art in spray computations in several important ways. 

15. SUBJECT TERMS 


Spray combustion modeling; Atomization; Sprays; Reacting flow modeling; Flash and superheat vaporization; CFD 


16. SECURITY CLASSIFICATION OF: 


17. LIMITATION OF 

18. NUMBER 

19a. NAME OF RESPONSIBLE PERSON 




ABSTRACT 

OF 

STI Help Desk (email:help@sti. nasa.gov) 

a. REPORT 

u 

b. ABSTRACT 

u 

c. THIS 
PAGE 

u 

uu 

PAGES 

62 

19b. TELEPHONE NUMBER (include area code) 

301-621-0390 


Standard Form 298 (Rev. 8-98) 
Prescribed by ANSI Std. Z39-18 




































